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ABSTRACT 


An analytical technique is developed to solve nonlinear three-dimensional, 
transverse and axial combustion instability problems associated with liquid- 
propellant rocket motors. The Method of Weighted Residuals is used to deter- 
’v. mine the nonlinear stability characteristics of a cylindrical combustor with 

uniform injection of propellants at one end and a conventional DeLaval nozzle 
•* at the other end. Crocco’s pressure sensitive time-lag model is used to des- 

cribe the unsteady combustion process. The developed model predicts the tran- 
sient behavior and nonlinear wave shapes as well as limit-cycle amplitudes 
and frequencies typical of unstable motor operation. The limit-cycle ampli- 
tude increases with increasing sensitivity of the combustion process to pres- 
sure oscillations. For transverse instabilities, calculated pressure wave- 
forms exhibit sharp peaks and shallow minima, and the frequency of oscillation 
is within a few percent of the pure acoustic mode frequency. For axial in- 
stabilities, the theory predicts a steep-fronted wave moving back and forth 
along the combustor. 
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SUMMARY 


An approximate analytical technique has been developed for the solution 
of nonlinear three-dimensional, transverse and axial combustion instability 
problems that are frequently observed in liquid-propellant rocket motors. 

This theory is an extension and generalization of previous analyses, which 
could analyze either transverse or axial instabilities in liquid combustors 
with quasi-steady nozzles, to the practical situations of three-dimensional 
instabilities in combustors with conventional DeLaval nozzles. Unlike the 
quasi-steady nozzle, the presence of a conventional nozzle imposes restric- 
tions upon the behavior of both the amplitudes and phases of the oscillations 
at the nozzle entrance plane. The Method of Weighted Residuals is used to 
determine the nonlinear stability characteristics of a cylindrical combustor 
with uniform injection of propellants at one end and a conventional nozzle 
at the other end. Crocco's pressure sensitive time-lag model is used to 
describe the unsteady combustion process. The developed model can predict 
the transient behavior and nonlinear wave shapes as well as limit-cycle ampli- 
tudes and frequencies typical of unstable motor operation. These results 
establish the relationship that exists between the resulting instability (i.e., 
waveform, final amplitude and final frequency) , the combustion parameters 
(i.e., interaction index, n, and time-lag, t) , and the chamber Mach number 
and length-to-diameter ratio. Results indicate that the limit-cycle ampli- 
tude increases with increasing sensitivity of the combustion process to pres- 
sure oscillations. For transverse instabilities, calculated pressure waveforms 
exhibit sharp peaks and shallow minima, and the frequency of oscillation is 
always within a few percent of the frequency of one of the chamber's acoustic 
modes. For axial instabilities, the theory predicts the presence of a steep- 
fronted wave moving back and forth along the combustor. In both cases calcu- 
lations of pressure and velocity perturbations at the nozzle entrance plane 
show that the approximation to the nozzle boundary condition is very good. 

The theory described in this report represents the final stage in the develop- 
ment of a unified nonlinear theory for the solution of general three-dimen- 
sional, transverse and axial combustion instability problems. 


INTRODUCTION 


Observation of the behavior of unstable rocket motors indicates that 
combustion instability can be divided into two categories; that is, linear 
and nonlinear instabilities. Linear instabilities are spontaneous in nature, 
and they are usually an outgrowth of the random combustion and flow fluctua- 
tions present in the system. On the other hand, nonlinearly unstable motors 
require the introduction of a finite amplitude disturbance to produce (or 
trigger) combustion instability. In either case the instability, after a 
transient period, reaches a limiting maximum amplitude (i.e., limit-cycle 
amplitude) at which it oscillates with a given frequency that is usually 
close to the frequency of one of the chamber's acoustic modes. Pressure 
measurements taken during test firings of unstable motors indicate that the 
limit-cycle waveforms of transverse instabilities are non -sinusoidal, that is, 
they exhibit sharp peaks and flattened minima."*" On the other hand, experi- 
mental observations of axial instabilities indicate the presence of shock- 
like steep-fronted waves in the chamber.^ These results indicate that non- 
linearities need to be considered in the theoretical treatment of combustion 

instability. 

Any analytical treatment of combustion instability should be capable of 
solving no nli near multi -dimensional combustion instability problems without 
exceeding memory core limitations of current computers and without requiring 
excessive computation time. To be of practical use, such a solution technique 
should be conceptually simple and easily adaptable for use by industry. This 
report describes the development and use of such a numerical solution tech- 
nique . 

Work on this problem has been in progress during the past several years , 
and due to its complexity, the problem had to be tackled in stages. In ear- 
lier investigations by these authors theories describing the nonlinear beha- 
vior of longitudinal^’^ and transverse^ instabilities in liquid combustors 
with quasi-steady nozzles were developed. These theories, which were based 
upon the application of the Method of Weighted Residuals (MWR) , successfully 


predicted the transient behavior, nonlinear -waveforms, and limit-cycle ampli- 
tudes of longitudinal and transverse instabilities in unstable liquid rockets. 
This report is concerned with the development of a generalized nonlinear 
theory that will be capable of analyzing three-dimensional, transverse and 
axial instabilities in the more practical situations where the combustors 
are attached to conventional nozzles. Obviously, this generalized theory 
will encompass the above-mentioned investigations as special cases. Contrary 
to the quasi-steady nozzle case, the presence of a conventional nozzle imposes 
both a mp litude and phase boundary conditions that must be satisfied by the 
solutions of the problem at the nozzle entrance plane. The generalized the- 
ory presented herein also provides a better description of the -unsteady flow 
field in the vicinity of the nozzle entrance plane. 

The application of the theory presented herein will be demonstrated 
by considering the nonlinear stability of a liquid-propellant rocket combus- 
tor with uniform injection of propellants at one end and a conventional noz- 

7 

zle at the other end. Crocco's pressure sensitive time lag model is used to 
describe the unsteady combustion process. In the sections to follow, the de- 
velopment of the wave equation for the analysis of nonlinear combustion in- 
stability in liquid rockets will be briefly described, the solution of this 
nonlinear wave equation will be outlined, and typical results will be present- 
ed and discussed. User's Manuals and program listings for the computer pro- 
grams used to solve these problems are included as appendices to this report . 


SYMBOLS 



*J‘i' I W t > 

time -dependent amplitudes in series given by Eq. 

(6) 

V 

A p (t) 

time -dependent amplitudes in series given by Eq. 

(9) 

r 

B(i) 

boundary residual 



b. 

Lmn 

complex axial acoustic eigenvalue 



* 

c 

velocity of sound, ft/sec 



3 


C 0’ C l’ C 2’ C 3 

D i 3 V V \ 

e(d 


t, m 


n 


P 


r 


coefficients of linear terms in Eqs.(l2) 
coefficients of nonlinear terms in Eqs . (12) 
residual of Eq. (10) 
imaginary unit, /- 1 

Bessel function of the first kind, order m 

axial and tangential mode numbers, respectively 

pressure interaction index 

* , * *2 

dimensionless pressure, yp /P 0 C Q 

* * 

dimensionless radial coordinate, r /R c 


c 

Rp (r) 

s 

ran 

t 


V 

•4 

Y 


chamber radius, ft 

radial acoustic eigenfunction in Eq. (9) 

dimensionless transverse mode frequency 

* 


dimensionless time, 


(R /c ) 
v c' o 7 


* . * 


dimensionless axial Velocity, u /c 


* / 

dimensionless velocity vector, V /c 


unsteady combustion mass source 


complex nozzle admittance 


* 

o 


4 


z 


Z. (z), Z (z) 
ton ’ p 


e 


e 

©(e) 

ir 

P 

T 


$ 


Subscripts: 

e 


n 


r, t, z, 0 


r, i 


o 


Superscripts : 


* , * 


dimensionless axial coordinate, z /r 


axial acoustic eigenfunctions 


ratio of specific heats 


ordering parameter 


azimuthal coordinate 


tangential acoustic eigenfunction in Eq. (9) 


•Wr , 

dimensionless density, p /p 


* 


dimensionless pressure sensitive time lag, — - — — 

(R /c ) 
v c' o 

velocity potential 


evaluated at the nozzle entrance 
radial mode number 

partial differentiation with respect to r, t, z, or Q 
respectively 

real and imaginary parts of a complex quantity, respec- 
tively 

stagnation quantity 


/ 


perturbation quantity, differentiation with respect to 
argument 


5 


steady state quantity 


* 




dimensional quantity, complex conjugate 
approximate solution 
ANALYSIS 


Development of the Wave Equation 

To keep the problem as simple as possible, yet still physically mean- 
ingful, the following assumptions are made. The gas phase in the combustor 
is assumed to consist of a single constituent which is thermally and caloric - 
ally perfect. Transport phenomena, such as diffusion, viscosity, and heat 
conduction are neglected. The momentum interchange between the liquid and 
gas phases is neglected (see Appendix A for a discussion of this assumption), 
and the specific stagnation enthalpy of the unburned propellant is assumed 

constant throughout the chamber . The presence of burning propellant drops 
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is represented by a distribution of unsteady mass sources and it is also 
assumed that the Mach number of the combustor's mean flow is small and that 

the waves have moderate amplitudes . 

As a result of the last two assumptions, the governing conservation 
equations ma y be combined and the unsteady flow in the combustor can be de- 
scribed by a single nonlinear wave equation. The derivation of this equation 
appears in Refs. 8 and 9, where it was assumed that each perturbation quanti- 
ty and the mean flow Mach number were of 0(e)> where e is an ordering para- 
meter that is a measure of the wave amplitude. After neglecting all terms of 
0(e 3 ) or higher and combining equations , one obtains the following nonlinear 
partial differential equation that describes the behavior of the velocity po- 
tential, $, of the combustor disturbance: 

V 2 § - * tt = + v(v.V)§ t + 2v$.V$ t + (y “ i) $ t v * + t 1 ) 

Equation (l) is the desired wave equation, and it is similar to the inhomo- 
geneous wave equation solved by Maslen and Moore 10 in a related study on non- 
linear acoustics. This equation accounts for the following effects: (l) the 
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effect of a steady state flow on the wave motion (viz., the first two terms 
on the right-hand side), (2) the coupling between the gas dynamical oscilla- 
tions and the unsteady combustion process (viz., the last term on the right- 
hand side), and (3) the second order nonlinearities of the gas dynamical 
processes (viz., the third and fourth terms on the right-hand side). 

In addition to satisfying Eq. (l), the desired solutions must satisfy 
rigid wall boundary conditions at the injector end of the chamber and at the 
chamber walls, while a nozzle admittance condition must be satisfied at the 
nozzle entrance. These boundary conditions are given (in a cylindrical coor- 
dinate system) by: 

$ = 0 at r = 1 

r 

$ = 0 at z = 0 

z 

B($) = yY^^ = ^ 2 = 

The nozzle admittance, Y, is a complex number defined by 

Y = Y + iY = (u7p') z = (3) 

r 1 e 

where u* is the dimensionless axial velocity perturbation and p is the di- 
mensionless pressure perturbation. 

It should be pointed out that due to the absence of an appropriate 
nonlinear nozzle admittance boundary condition, the solutions of the problem 
are required to satisfy a linear nozzle admittance. Although inconsistent 
with the nonlinear wave equation, the linear nozzle admittance condition is 
used herein with the hope that the solution techniques developed herein will 
also be applicable when nonlinear nozzle admittance conditions become avail- 
able. Also, the relative importance of nozzle nonlinearities is not known 
at the moment and it is quite possible that the linear nozzle boundary condi- 
tion used herein adequately describes the flow conditions at the nozzle en- 
trance . 

The unsteady combustion process is represented by mass sources distri- 
buted throughout the volume of the chamber, and the response of the mass 
sources to pressure oscillations is assumed to be described by Crocco's pres- 
sure sensitive time-lag hypothesis.^ The mass source perturbation, W m , 

5 8 

is then given by: ’ 
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(4) 


W' = -vn 
m 1 


du r 
dz L 


§ f (r,e,z,t) - $ t (r,9,z,t - f) 


where n is the pressure "interaction index" that describes the sensitivity of 
the combustion process to pressure oscillations, and j, commonly referred to 
as the sensitive time-lag, is the part of the total combustion time-lag during 
which the combustion process is sensitive to pressure oscillations . The un- 
steady combustion response described by Eq. (4) is linear and the com me nts 
made above regarding the use of a linear nozzle admittance boundary condition 
are also applicable to this case. 

Substituting Eq. (4) into Eq. (l) and expressing the resulting equation 
in a cylindrical coordinate system yields the following wave equation: 


(5 4 " — — $ -|- * ■■ $ + 

®rr r *r 2 ^99 
r 


zz 


tt 


- 2 $ $ . - — $ - 2 $ $ 
r rt 2 0 9t z zt 

r 

- (V - 1)4 (* + 7 $ + -4 + $ ) 

t rr r r r 2 99 zz 
. du 

“ 2u *zt " Y *t dl 


+ yn dz “ 4 t (r,9,z,t - t)] = 


(5) 


The combustor and nozzle geometries considered in this study, as well as the 
cylindrical coordinate system used in writing Eq. (5), are shown in Fig. 1. 

Method of Solution 

Since Eq. (5) has no known closed-form mathematical solution, it is 
necessary to resort to the use of either exact numerical solution techniques 
or approximate analytical techniques. For multi -dimensional problems, the 
exact numerical solution techniques generally exceed the computer storage 
capacities, therefore an approximate solution technique is used herein. 

The experience of previous investigators in the fields of structural stabili- 
ty and aeroelasticity indicates that an approximate solution technique known 

11 12 

as the Method of Weighted Residuals 5 may be effective in the solution 
of this nonlinear wave equation. 

In order to employ the Method of Weighted Residuals in the solution 
on Eq. (5), it is first necessary to express the velocity potential, I, as an 
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Figure 1. Combustor Configuration and Coordinate System. 


approximating series expansion, §. The question naturally arises as to what 
form of series expansion should he used. Inasmuch as the experimentally ob- 
served pressure oscillations during combustion instability usually resemble 
the natural acoustic modes of the chamber, the velocity potential, $, is ex- 
panded in terms of the natural acoustic modes of the chamber with unknown 
time -dependent amplitudes. 

In previous analyses' 5 ’ of related problems the approximate solutions 
were expressed in terms of the acoustic modes for a cylindrical chamber with 
solid wall boundary conditions at both the injector and the nozzle ends. 
Consequently, the approximation of the flow conditions at the nozzle entrance 
was poor. In the present analysis a better approximation to the flow at the 
nozzle entrance is obtained by expanding the velocity potential in terms of 
the acoustic eigenfunctions for a chamber with a solid wall boundary condition 
at the injector end and a nozzle admittance condition at the other end. This 
removes both the two -dimensionality and the quasi -steady nozzle restrictions 
imposed upon the previous investigations . 

The velocity potential, § , is therefore approximated by the following 
series expansion: 

% = Y V Y {a (t) sin me + B, (t) cos me} Z (z) J (S r) ( 6 ) 

/././. L J-mn ton J ton m mn 

t m n 

where the A’s and B’s are unknown complex functions of time, and the Z's are 
the complex axial acoustic eigenfunctions. The complex form of the axial 
acoustic eigenfunctions is given by 


Z„ (z) = cosh(ib, z) 
tnur ton 

where the b. are the axial acoustic eigenvalues which must satisfy the 
ton 

following transcendental equation: 


(7) 


b L sin2 <’We ) 


+ y^Y 2 (S 2 + b 2 ) cos 2 (b 
T v mn ton v 


, z ) = 0 
ton e 


(8) 


Equations (7) and ( 8 ) are obtained by linearizing Eq. ( 5 ) and solving the re- 
sulting equation for the case of no mean flow or combustion (i.e., the acous- 
tic case) subject to the boundary conditions specified in Eq. (2) . Each 
term in the above expansion exactly satisfies the solid wall boundary condi- 
tions at the injector end (i.e., at z = 0 ) and at the chamber wall (i.e., at 
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r - 2 .) ; however, due to the unknown time dependence of Eq. (6) the nozzle 
admittance condition imposed at z = z g is not exactly satisfied by the indi- 
vidual terms . Including both the sin me and cos m 9 terms in the expansion for 
allows for the possibility of either spinning or standing wave solutions . 

In order to simplify the algebra involved in the application of the 
Method of Weighted Residuals, the development of the associated computer 
program, and the presentation of the results; the expansion of the velocity 
potential is written as a single summation as follows: 

N 

i = Y A (t)Z (z)© (e)R (r) (9) 

L P P p P 

p=l 

where the A ’ s are the unknown time -dependent amplitudes. In order to use 

p 

Eq. (9) a correspondence must be established between the index, p, in Eq. (9) 
and the mode-numbers l, m, and n in Eq. ( 6 ). Such a correspondence is given 
in Table 1 for a three mode series consisting of the spinning first tangential 
(IT) mode (t = 0 , m = 1 5 n = l), the spinning second tangential (2T) mode 
(1=0, m = 2, n = l), and the first radial (lR) mode (<& = 0, m = 0, n = l) . 


Table 1 

Correspondence Between Eq. ( 6 ) and (9) for a Three-Mode Series 


p 

Mode 

•L(p) 

m(p) 

n(p) 

A 

P 

0 

J2 

1 

IT 

0 

1 

1 

A 0 u(t) 

sin 9 

2 

IT 

0 

1 

1 

B 011 (t) 

cos 9 

3 

2T 

0 

2 

1 

A (t) 
021 

sin 29 

k 

2T 

0 

2 

1 

B 021 (t) 

cos 29 

5 

IE 

0 

0 

1 

B oof ^ 

1 


Before proceeding with the analysis, the wave equation (i.e., Eq. (l)) 
must be modified for use with the assumed complex solution given by Equation 
(9) • This modification is necessary because only the real part of the assum- 
ed solution is physically meaningful. It can easily be shown that if $ = cp 
+ iY is a solution to Eq. (l), the real part, 9 , is not a solution to Eq. (l) . 
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This failure of <p to satisfy Eq. (l) is due to the presence of the nonlinear 
terms in this equation. It can also be shown, however, that a modified wave 
equation can be constructed for which the real part of its solution satisfies 
the original wave equation (i.e., Eq. (l)). This modified wave equation is 
given by: 

E(t) = V 2 $ - - 2V . v« t - Y(V.V) * t - W m 

- ■ - 1 ' — j^v4.V« t + V§*v£*] - [v$.V^ + V$*.V$ t ] 

- ■ f — {(1 - i) [$ t V 2 S + $*V 2 $*] 

+ (l + i) ^$^v 2 §* + $*V 2 |J 1 = 0 (10) 

where $* is the complex conjugate of $. The derivation of this equation is 

discussed in Appendix B. Thus, the Method of Weighted Residuals will be used 

to obtain approximate solutions to Eq. (10) (i.e.,?=cp+ iY) from which the 

real part, will be taken as the approximate solution of Eq. (l) . 

In order to obtain a solution, the unknown time -dependent mode -amplitudes 

(i.e., A (t)) are determined by the following mathematical procedure. The as- 
P 

sumed series expansion, (i.e., Eq. (9)) is substituted into the wave equation 
(i.e., Eq. (10)) to form the equation residual, E('l) . Similarly, substituting 
the series expansion into the nozzle boundary condition (i.e., the last of Eq. 
(2)) yields the boundary residual, B($) . In the event that these residuals 
are both identically zero, the solution is an exact solution. The residuals 
E($) and B($) represent the errors incurred by using the approximate solution, 

According to the modified version of the Method of Weighted Residuals, 
developed by the authors in Refs. 5 and 8, the residuals E(l>) and B($) must 
satisfy the following orthogonality conditions: 

Z e 2 tt 1 

J J J E(!)z*(z)0^(e)R^(r)rdrdedz 
0 0 0 

2n 1 

- f f B($)Z*(z )® . (e)R . (r)rdrd0 =0 
J J 3 ® 3 3 

0 0 


( 11 ) 


where in the present study the complex conjugate of the axial eigenfunction, 

Z*, is used in the weighting functions. The chosen weighting functions must 
correspond to the terms that appear in the assumed series solution; that is, 

Eq. (9). 

Evaluating the spatial integrals in Eq. (ll) yields the following system 
of E complex ordinary differential equations to be solved for the unknown com- 
plex amplitude functions, A (t) : 


( d A r 

^ {C Q (j,p) \ + c 1 (j ,p)A p (t) + [C 2 (j,p) 

p=l 



dA 

dt 


d[A (t - t)] 
+ nC 3 (j,p) 


N N 


dA, 


dt 


} + I I A p “at + V 3 ’ p >' l)A : 

p=l q=l 


p dt 


* ^Aq * 

+ Vd.P.OA It + D 4 (3 ’ p>q)A p “dt ; = 0 


( 12 ) 


j — 1 , 2 , ... N 


The coefficients appearing in the above equations are determined by evaluating 
the various integrals of hyperbolic, trigonometric, and Bessel functions that 
arise from the spatial integrations indicated in Eq. (ll) . A user's manual for 
the computer program C0EFFS3D used to calculate these coefficients is given in 
Appendix C . 

The time-dependent behavior of an engine following the introduction of 
a disturbance is determined by specifying the form of the initial disturbance 
and then following the subsequent behavior of the individual modes by numerical 
ly integrating Eqs. (12). Once the time -dependence of the individual modes is 
known, the velocity potential, §, is calculated from Eq. (9)* -*-he pressure 

perturbation at any location within the chamber is related to the real part of 
¥ (i.e., cp) by the following second-order momentum equation (see Refs. 5 and 8) 


p' = -YtqJt + u(z)9 z + + - 2 <Pe 


rJ2. ~ 2 , 


!~2- 


(13) 
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A user's i narm a. 1 for the computer program, LCYC3D, which obtains numerical so- 
lutions of Eqs. (12) and (13) is given in Appendix D. 

In summary, the theory presented in this section represents a two-stage 
simplification of the original problem. In the first stage the problem has been 
reduced to the solution of a single nonlinear, partial differential equation 
(i.e., Eq. (l)). In the second stage the solution was expanded in a series of 
acoustic modes with time -dependent coefficients and the Method of Weighted 
Residuals was used to replace the solution of the nonlinear partial differential 
equation with the solution of a system of nonlinear, ordinary differential 
equations (i.e., Eq. (12)). Typical numerical solutions of these equations will 
be presented and discussed in the following section. 

RESULTS AND DISCUSSION 

The generalized three-dimensional theory introduced in the previous 
section has been used to obtain both linear and nonlinear data for pure trans- 
verse modes and pure longitudinal modes for rocket motors with conventional 
nozzles. Nonlinear data for the first tangential (IT) mode and the first lon- 
gitudinal (li) mode has also been obtained for combustors with quasi-steady 
nozzles for comparison with the results of the previous two-dimensional theories. 
Linear Solutions 

Before proceeding with the nonlinear analysis, it was desired to obtain 
numerical solutions of the linearized equations (i.e., Eqs. (12) with D^ = D^ = 

= D^ = 0) in order to determine how closely the approximate solutions satis- 
fied the nozzle boundary condition. The linear solution is also needed for com- 
parison with the corresponding nonlinear results. The linear solutions were 
obtained by assuming a one-mode series expansion consisting only of the mode 
under consideration. Due to the presence of the retarded variables (i.e., 

d[A (t - T)]/dt) in Eqs. (12), it is necessary to specify the initial amplA- 
P 

tudes over the interval -t £ t £ 0. In this study the initial values were 
chosen such that the nozzle boundary condition was exactly satisfied during 
this initial time period. Solutions were obtained for values of n and t on 
the neutral stability limit (see Appendix E for the determination of neutral 
stability limits) for various conventional nozzle configurations. The nozzle 
admittance was expressed in the form, Y = Ae lc P, where A is the amplitude fac- 
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tor and cp is the phase shift. The pressure perturbation, p' , ana the axial 

velocity perturbation, u' , at the nozzle entrance were calculated numerically 

for several values of the nozzle phase shift, cp. These calculated values were 

then used to compute the ratios (u / /p , )„_„ > which were then compared with the 

e 

snecified nozzle admittance values. These results are shown in Tables (2) and 
( 3 ) where A and cp n are the computed values of the amplitude factor and phase 
shift, respectively. These results show that the approximation to the nozzle 
boundary condition is very good for both the IT and 1L modes; that is, the 
maximum error in the amplitude ratio is about 5 °jo and the maximum error in 
phase is approximately 0.5 degree. These results are in contrast with previous 
theoretical investigations where the representation of the unsteady flow con- 
ditions in the vicinity of the nozzle entrance was very poor. 

Table 2. IT Mode Linear Solutions (Numerical). 



A = 0.02 


Error at 

Nozzle 

cp 

(Degrees) 

T 

a 

A - A 
n 

A 

9 n - 9 

(Degrees) 

0 

1.2 

0 . 664l6 

-.029 

0.4 

0 

1.7 

0.55001 

.003 

0.4 

0 

2.2 

0 . 64710 

.034 

0.4 

45 

1.2 

0.66137 

-.031 

0.5 

45 

1.7 

0 . 54490 

.001 

0.5 

45 

2.2 

0.63665 

.032 

0.5 

90 

1.2 

0.62507 

-.031 

0.3 

90 

1.7 

0.51252 

-.001 

0.3 

90 

2.2 

0.59758 

.028 

0.3 

135 

1.2 

0.57746 

-.031 

-0.1 

135 

1.7 

0.47274 

-.004 

-0.1 

135 

2.2 

0.55353 

.023 

-0 . 1 

180 

1.2 

0.54825 

-.030 

-0.4 

180 

1.7 

0.45003 

-.004 

-0.4 

180 

2.2 

0.53121 

.022 

-0.4 

225 

1.2 

0.55357 

-.030 

-0.4 

225 

1.7 

0.45677 

-.003 

-0.4 

225 

2.2 

0.54292 

.024 

- 0.5 
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Nonlinear Solutions 

Nonlinear solutions have been computed for both the IT mode and the LL 
mode. For the IT mode calculations a three mode series expansion consisting 
of the IT, 2T (second tangential), and 1R (first radial) modes was used. These 

are the sa me modes that were included in the series expansion used in the pre- 

5 6 

vious two-dimensional transverse instability studies. ’ In these studies it 
was shown that convergence was obtained with this three mode series; that is, 
the addition of higher transverse modes (i.e., 3T, 4T, etc.) to the basic series 
had little effect on the solution. The 1L mode computations were made using a 
series consisting of the first five longitudinal modes (i.e., 1L, 2L, 3L, 4h, 
and 5L) . It has been shown by Lores and Zinn ’ that convergence is obtained 
with this five -mode series. 

Transverse Mode Solutions . Nonlinear solutions have been computed for 
rocket combustors with quasi-steady nozzles (i.e., real admittances) and also 
for nozzles with complex admittances. The quasi-steady nozzle solutions were 
generated for comparison with the results of the previous two-dimensional the- 

5 13 

ory. For this case the nozzle admittance is given by: 



Y i = 0 (14) 

For nozzles with complex admittances the admittance was expressed in the form, 
Y = Ae lcp . For both cases limit -cycle amplitudes and waveforms have been com- 
puted for both standing and spinning first tangential instability. This re- 
quired three series terms to describe standing instability and five series 
terms to describe spinning instability. Typical computation times on a Uni- 
vac 1108 computer to reach a limit -cycle were one minute for a standing wave 
and two minutes for a spinning wave . 

Wall pressure waveforms (r = l) were computed at the injector face 
(z = 0) and at the nozzle entrance (z = z g ) for three azimuthal locations, 

9 =0°, 0 =45°, and 0 = 90°. The initial conditions for standing waves were 
chosen such that a pressure anti -node occurred at 9=0. Injector pressure 
waveforms for both standing and spinning instability are shown in Fig. 2 for 
combustors with quasi-steady nozzles. These waveforms exhibit sharp peaks 


17 


Dimensionless Time, t 


(b) Spinning Wave: n = O.58, t = 1.7, Y = 1.2, u = 0.2, z = 1 

e ’ e 


Figure 2. Nonlinear Pressure Waveforms for the IT Mode 








and shallow minima; they are nearly identical in shape to those calculated 

5 6 

using the previous two-dimensional theory. Comparison of injector and 

nozzle pressure waveforms (9 = 0°) shows that there is very little variation 
in pressure with axial position. These waveforms are in qualitative agree- 
ment with the results of pressure measurements taken during test, firings of 
unstable rocket motors. 

To check the accuracy of the approximation of the nozzle boundary 
condition, wall pressure and axial velocity waveforms were calculated at the 
nozzle entrance. The error at the nozzle boundary (z = z^) is shown for non- 
linear standing and spinning IT mode instabilities in Fig. 3* Here the axial 
velocity perturbation, u', and the product of the quasi-steady nozzle admit- 
tance and the pressure perturbation, Y^p 7 are plotted as a function of time. 

The latter quantity is the axial velocity perturbation that would be obtained 
at the nozzle entrance if the nozzle boundary condition were exactly satis- 
fied (i.e., the nozzle admittance condition requires that u 7 = Y^p 7 at z = z g ) • 
Most of the discrepancy between the two curves is due to a slight phase shift 
Pg^ w g 0 fi pressure and velocity and the second harmonic distortion of the pres- 
sure waveform resulting from the nonlinearities of the system. The nozzle 
boundary condition is satisfied in an average sense, however, for the ratio 
of the velocity amplitude (peak-to-peak) to pressure amplitude (peak-to-peak) 
is very close to the required value, Y^. 

In another study, limit-cycle amplitudes were calculated as a function 
of n and t for standing IT mode instability. Values of n in the linearly un- 
stable region were chosen for below resonant (t = l>9)s resonant (t = 1.706), 
and above resonant (t — 1.5) conditions. The resulting amplitudes are compared 
with those obtained with the two-dimensional theory in Fig. 4. This figure 
shows that the three-dimensional theory predicts a slightly higher limit-cycle 
amplitude than the two-dimensional theory for chambers with quasi-steady nozzles. 

Figure 4 also shows that the three-dimensional theory, like the previous 
two-dimensional one, cannot predict triggering of IT mode instability by the 
introduction of finite amplitude disturbances. This result was expected since 
it was shown in Refs. 6 and 8 that the second order (i.e., 0( e )) theory can 
predict triggering only for pure radial modes (m = 0, n = 1, 2 ...). Such 
triggering limits for the 1R mode are discussed in Ref. 9* It has also been 
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Figure 4. Limit-Cycle Amplitudes for the IT Mode 





shown, however, that triggering of H mode instability can be described when 
the 0(* :S ) terms are retained in the analysis. ’ The third order theory 
given in Refs . 8 and 14 is limited to a single mode in the approximating 
series expansions. A more general multi -mode third-order theory is now under 
development and the results will be presented in a future publication. It 
is also suspected that nonlinear unsteady combustion effects (not included in 
the present analysis) may play an important role in the triggering phenomenon. 

For nozzles with complex admittances a study was conducted to determine 
the effect of the nozzle phase shift, cp, upon the limit-cycle amplitudes and 
waveforms for both standing and spinning IT mode instability. The effect of 
nozzle phase shift on the nonlinear pressure and velocity waveforms at the 
nozzle entrance plane is shown in Fig. 5 for spinning waves. This figure 
shows that, while cp has little or no effect on the pressure waveforms, the 
phase and shape of the velocity waveforms is strongly dependent on cp. The 
effect of cp on the limit-cycle amplitude for standing IT mode instability is 
shown in Fig. 6. For a given value of n and t (in the linearly unstable re- 
gion for the IT mode), Fig 6 shows a sinusoidal variation of limit-cycle am- 
plitude with cp having a maximum amplitude at about cp = 200° and a minimum am- 
plitude at about cp = 20°. In this connection, it should be pointed out that 
according to linear results nozzle damping is a maximum at cp = 0° and a mini- 
mum at cp = 180°; thus the observed shifts must be due to nonlinearities. 

In order to determine how well the solutions approximate the nozzle 
boundary condition, the amplitude ratio and phase shift between pressure and 
velocity at the nozzle entrance have been calculated from the nonlinear solu- 
tions and have been compared with the specified nozzle admittance condition. 
Since the waveforms are non -sinusoidal, an approximate amplitude ratio, A c , 
was calculated by taking the ratio of peak -to -peak velocity amplitude to 
peak-to-peak pressure amplitude. The approximate phase shift, cp c was calcu- 
lated from the following formula: 

cp c = p- ~ ty ] x 360 (15) 

where t is the average of an ascending zero-crossing and the following de- 

ir 

scending zero-crossing for the pressure perturbation, t is a similar average 
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Nozzle Axial Velocity 





Wall Pressure Amplitude (pk-pk) 






for the velocity perturbation, and T is the period of oscillation. The re- 
sults of this study are shown in Fig. 7 for both standing and spinning waves. 
For standing waves the calculated amplitude ratios are seen to be consistent- 
ly higher than required by the nozzle admittance condition (dashed line) , 
while for spinning waves the calculated amplitude ratios are lower than re- 
quired. For both standing and spinning waves the calculated phase shifts are 
in excellent agreement with the imposed phase shifts. This study shows that 
the three-dimensional theory provides a good approximation to the nozzle 
boundary condition for the IT mode, considering that the nonlinear solutions 
are being forced to satisfy a linear boundary condition. 

Longitudinal Mode Solutions . Letting m and n equal zero in Eq. (6) 
and using a series consisting of the first five longitudinal modes (i.e., t = 
1, 2, ... 5)3 limit-cycle solutions were calculated for quasi -steady nozzles 
as well as for nozzles with complex admittances. The longitudinal mode so- 
lutions required somewhat longer computation times than the transverse mode 
solutions; the time required to reach a limit cycle was from three to four 
minutes on the Uni vac 1108 computer. 

Longitudinal mode solutions for chambers with quasi-steady nozzles 

3 4 

were compared with the solutions previously obtained by Lores and Zinn ’ 

using a one -dimensional theory. Pressure waveforms at the injector face 

are compared for both resonant and off -resonant conditions in Fig. 8 which 

shows excellent agreement between the two theories. Pressure and velocity 

waveforms at the nozzle entrance as well as injector face pressure waveforms 

are shown in Fig. 9 for quasi-steady nozzles, while Fig. 10 shows waveforms 

at the nozzle entrance for nozzles with complex admittance (cp = 45° and 

q> = 90 °) • In each case the results indicate the presence of a steep-fronted 

pressure wave moving back and forth in the chamber. This behavior is in 

2 

agreement with experimental observations of axial instabilities. The 
relation between pressure and velocity waveforms at the nozzle entrance is 
a fairly good approximation to the nozzle admittance condition (see Figs. 9 
and 10) in spite of the highly nonlinear waveforms. The results of this 
investigation indicate that the three-dimensional nonlinear theory is appli- 
cable to longitudinal instabilities as well as transverse instabilities. The 
theory can also be uied to investigate the nonlinear behavior of combined 
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Figure 10. Longitudinal Mode Waveforms for Nozzles 
with Complex Admittances. 
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longitudinal-trans verse instabilities, although no results for instabilities 
of this type are presented in this report. 

CONCLUDING REMARKS 

A general three-dimensional second-order nonlinear theory has been 
developed for predicting the linear and nonlinear behavior of combustion in- 
stability in liquid-propellant rocket combustors. This theory contains previous 
analyses of transverse and longitudinal instabilities as special cases. Further- 
more it extends the previous analyses which were applicable only to combustors 
with quasi -steady nozzles, to the more practical cases of combustors with con- 
ventional DeLaval nozzles. The present theory can be used to predict the sta- 
bility characteristics of longitudinal, transverse and combined longitudinal- 
transverse modes for various liquid-propellant rocket motor designs. 

Results obtained for combustors with quasi-steady nozzles are in excel- 
lent agreement with the predictions of previous theories for both transverse 
and longitudinal instabilities. For combustors with conventional nozzles the 
limit-cycle amplitude varies sinusoidally with nozzle phase shift, cp, having a 
maximum value at cp = 200° and a minimum value at cp = 20°. The nozzle phase 
shift has a strong effect on the axial velocity waveforms at the nozzle entrance 
while having only a minor influence on the nonlinear pressure waveforms . In 
both cases, the nonlinear theory developed in this paper provides a good approx- 
imation to the unsteady flow conditions at the nozzle entrance plane. This is 
in contrast to the previous theories which provided a relatively poor approxi- 
mation to the nozzle boundary condition. 

The results presented in this report establish the relationship that 
exists between the resulting instability (i.e., waveform, final amplitude, and 
final frequency), the combustion parameters (i.e., interaction index, n, and 
time-lag t) , and the chamber Mach number and length-to-diameter ratio. These 
results indicate that the limit-cycle amplitude increases with increasing sen- 
sitivity of the combustion process to pressure oscillations . For transverse 
instabilities, calculated pressure waveforms exhibit sharp peaks and shallow 
minima, and the frequency of oscillation is always within a few percent of the 
frequency of one of the chamber's acoustic modes. For axial instabilities, the 
theory predicts the presence of a steep-fronted wave moving back and forth 
along the combustor. In both cases the calculated pressure waveforms are in 
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good qualitative agreement with available experimental data. 



appendix a 


MOMENTUM INTERCHANGE BETWEEN LIQUID AND GAS PHASES 

that 7 reSUltS PreMnt6,i ln this report — obtained under the assumption 
t at the momentum interchange between the liquid droplets and the burned gases 

g agible. This assumption will now be relaxed for the special case of 
uni ormly distributed combustion, and it will be shown that this momentum 
interchange is an important stabilizing effect. 

Analysis 

The momentum equation for two-phase flow was derived in Ref 8 and is 
given by: b 


Sv 


i 3? + 1 ™ ] ♦ - -<v - ^Hc + , ) 


(A-l) 


l 7 i *" the f and 11 ^ id respectively. The term on th 

an si e of Eq. (A-l) represents a momentum source to the gas produced 

the t ™ ng ll,Uid dr ° PS ' ThlS " K '‘ KritU '" S ° UrCe oooeiete of two parts: (l 

the force necessary to accelerate the evolved gases from the droplet velocity 

to the gas velocity (i.e., the term -W (V - V )) and (2) the 

of the droplets (i.e., the term -C(V -£>) * —dynamic drag 

is necl° rde 7 0 deriTC " ^ e9Ua " ion ^ OT Telocity potential * it 

eary to make the following assumptions: (l) the drag term is negli- 

are negTigibl 7 aCCeleraU ° n ta ™> < 2 > ^ velocity fluctuations 

chamber. Neglecting theT “ Unif0r,,,ly dlstl ' ita ted throughout the 

titles gives the Z ’ ^ Neglecting third order quan 

following expression for the momentum source perturbation M- 


5' “ - ( r - l'>«m (A-2) 

This is simplified further by neglecting the liquid velocity perturbation, 
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introducing the velocity potential, and using the steady-state relation, 

¥ = du/dz, to obtain: 

m ' 


M' = - f VJ (A-3) 

Finally, the assumption of uniformly distributed combustion gives du/dz = 
constant which yields: 


M' 

-¥ 



(A-4) 


Perturbing the left-hand- side of Eq. (A-l), introducing the velocity potential, 
and combining with Eq. (A-4) gives: • 

l3t + V + 5 h + i 4 + t 7} ' 7 *-K] = 0 (A ' 5) 

which can be integrated to obtain: 

p' = -r [ $ t + u$ z + fz $ + I V$ ‘ V$ - \\ ] (A-6) 


Equation (A-6) is similar to Eq. (13), where the additional term (du/dz) $ 
arises from the droplet momentum source. Following the procedure outlined 
in Ref. 8, the momentum equation given by Eq. (A-6) is combined with the con- 
tinuity and energy equations to obtain the desired wave equation: 


V 2 $ 


$ tt ~ 2S$ zt 


+ Df *t + 


2V$-V$ t + (r - l)$ t 7 § + w m ' (A-7) 


Comparing Eq. (A-7) with Eq. (l) shows that the droplet momentum source 
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appears only in the second term on the right -hand- side of this equation, 
where the factor r in Eq. (l) becomes (r + l) in Eq. (A-7) . 

Applying the Method of Weighted Residuals to obtain approximate solutions 
to Eq. (A-7) yields a set of ordinary differential equations identical to Eqs. 
( 12 ) where the coefficient Cg(d,p) is now given by: 


z z g 

»(a.P> - { 2 / e 5(z)z p ^dz + (r + 1)J 0 fz p z*dz + nz p (z e )z*(z e )} x 

X \ © 0 .d 0 |" R R.rdr (A- 8 ) 

» „ P J J Q PJ 


Equation (A- 8 ) is readily obtained from Eq. (C-3) by replacing r in the second 
term by y + 1 . 

Linear Stability Limits 

Linear stability limits for the 1L mode were calculated by the method 
described in Appendix E for the following two cases: (l) the droplet momen- 

tum source was included in the analysis and ( 2 ) the droplet momentum source 
was neglected. The results were compared with the linear stability limit 
calculated by Mitchell 1 ^ on a plot of interaction index, n, versus stretched 
time-lag, p,, where p, = ot/tt (see Fig. A-l) . This figure shows excellent agree- 
ment between the results of Mitchell (solid curve) and the present theory 
(circle symbols) when the droplet momentum source is included. Neglecting the 
droplet momentum source shifts the stability curve to much lower values of n 
(dashed curve) , which indicates that the droplet momentum source is an impor- 
tant stabilizing effect. 

Nonlinear Solutions 

In the second-order analysis presented in this report, the droplet momen- 
tum source affects the nonlinear solutions primarily by increasing the linear 
stability of the system. This is readily shown in Fig. (A-2) where the limit- 
cycle amplitude is plotted as a function of the displacement, fin, above the 
neutral stability limit. By plotting the limit-cycle anplitudes in this manner 



Interaction Index 



Figure A-l. Effect of Droplet Momentum Source on Linear Stability- 
Limits for the 1L Mode. 
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the effect of the shift in the neutral stability curves is removed so that 
only the nonlinear effect of the momentum source is seen. Figure A-2 shows 
that, for equal displacements above the neutral stability limits, including 
the droplet momentum source results in a slightly smaller limit-cycle ampli- 
tude. This difference in limit-cycle amplitude is negligible for most practi- 
cal purposes. 

For combustors with uniformly distributed combustion it has been shown 
that the droplet momentum source is an important effect which is easily in- 
corporated into the present analysis. Consequently the computer programs 
based on this theory include the droplet momentum source as an optional feature 
(see Appendices C, D, and E) . 

For chambers with non-uniform combustion distributions, Eqs. (A- 6 ) and 
(A-7) are no longer applicable; however, the droplet momentum source can be 
taken into account in the following manner. Using the present theory with 
the droplet momentum source omitted, the neutral stability limit, n (t), is 
calculated and the limit-cycle amplitudes are determined as a function of 6 n 
as in Fig. A-2. In addition, the linear stability limit, n (t), is calculated 
using a linear theory which includes the droplet momentum source and is not 
restricted to uniformly distributed combustion (such as in Ref. ( 15 )). Assum- 
ing that the nonlinear effect of the droplet momentum source is also small for 
non-uniformly distributed combustion and using the values of 5 n and n^(j) cal- 
culated above, the desired plot of limit-cycle amplitude as a function of n 
is readily obtained. 


37 


APPENDIX B 


USE OF COMPLEX VARIABLES IN THE SOLUTION OF 
NONLINEAR DIFFERENTIAL EQUATIONS 

It is often convenient to use complex variables in the solution of the 
linear equations which arise in acoustics, combustion instability and related 
fields. In this case the solution is expressed in complex form, and the real » 
part represents the physically meaningful solution. However, care must be 
used when applying this technique in the solution of nonlinear equations. The 
difficulties that are encountered in applying the complex variable technique 
to nonlinear problems will be illustrated by analyzing the following simplified 
example. Consider the nonlinear wave equation given by: 

v 2 $ - * tt = ii t (B-l) 

A complex solution of Eq. (B-l) of the form § = q> + iy would be useful only 
if its real part, qj, satisfies Eq. (B-l), which would be the case if the equa- 
tion were linear. However, straightforward substitution of $ = m + iy into 
Eq. (B-l) and separating its real and imaginary parts yields the following 
equation for qj : 


v2 9 “ <Ptt " Wt " n t ( B ~ 2) 

indicating that the real part, qj, does not satisfy Eq. (B-l) because of the 
extra term, -yy , appearing on the right hand side. In order to eliminate 
this extra term, the form of the original differential equation (i.e., Eq. 
(B-l)) must be modified. 

Since Eq. (B-l) supposedly describes some physical phenomenon, and 
since only the real part of the complex solution is physically meaningful, 
then the nonlinear term should really be expressed as the product Re($) X 

^ -X- 

Re($^.) which is equivalent to ($f^ + + $ $^)/4. Substituting this 

expression into Eq. (B-l) yields: 
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2 -x- -x- X" 

V $ - $ tt = + $$ t + $ $ t + $ $ t ] 


(B-3) 


Substituting f = qj + iy into Eq. (B-3) and separating its real and imaginary 
parts yield: 


2 

V <p - cp tt - cpcp t 

(B-4) 

v 2 y - Y tt = 0 


which shows that the real part of the solution of Eq. (B-3) satisfies the de- 
sired equation (i.e., Eq. (B-l)) and the imaginary part satisfies a homogeneous 
linear wave equation. This technique was applied to the solution of nonlinear 
combustion instability problems (i.e., to Eq. (l)), and the resulting modified 
wave equation was solved using the Method of Weighted Residuals. Due to the 
approximate nature of the Method of Weighted Residuals, however, the resulting 
solution contained an error term which grew without limit. Consequently, the 
above procedure had to be modified in order to obtain satisfactory solutions 
of Eq. (l) using the Method of Weighted Residuals. 

An alternate technique is to modify Eq. (B-l) such that both the real 
and imaginary parts satisfy the original equation. This can be done by re- 
placing terms of the form with Re($)Re($ i _) + ilm($)lm($^) ; using the 
relations: 


Re($)Re($ t ) = ( § $ )( — — t ) = + #*§*J 


- * 




(B-5) 


in Eq. (B-l) gives: 

2 f“ X- "X" ^ y — 

V § - § tt = £[_(! - i)(f§ t + § $ t ) + (l + i)($$ t + $ § t )J (B-6) 

Substituting $ = cp + iy into Eq. (B-6) and separating into its real and imagi- 
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nary parts gives: 


2 

v <P " «Ptt ~ Wt 


- T tt = ff t 


(B-7) 


which shows that both <p and y satisfy Eq. (B-l). Applying this method to the 
solution of Eq. (l) yields the modified wave equation (i.e., Eq. (lO)) used in * 
the present investigation. 
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APPENDIX C 

PROGRAM COEFFS3D: A USER'S MANUAL 


Statement of the Problem 


Program C0EFFS3D calculates the coefficients of both the linear and non- 
linear terms which appear in Eqs . (12). These coefficients are required as 
input for Program LCYC3D (see Appendix D) which numerically integrates this sy- 
stem of equations. The coefficients that are required depend on the choice of 
terms to be included in the series solution for $ (see Eq. (9)), therefore 
this information must be provided as input to Program C0EFFS3D. The output 
of Program C0EFFS3D is either punched onto cards or stored on drum (FASTRAND) 
for input to Program LCYC3D. 


The coefficients to be calculated are functions of various integrals 
of hyperbolic, trigonometric, and Bessel functions and are given by the follow- 


ing expressions : 
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2tt 


(j,p) = f z Z* dz f © 0.d0 f 
« P J J P J J 
0 0 


R R . rdr 
P J 


(C-l) 
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C 1 (j,p) 


{sjj(p) f z Z* dz - f z"Z* dz + Z'(z )Z*(z )} f @ O.de f R R.rdr 
L ran J p 0 J p 3 P e' e'J J p j J p j 
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c 2 (3 ,p) 


{2 J i(z)Z'Z* a, , V J | ^ dz + Y*Z p < Z e > Z > e >} 
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“ PJ J P J 
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= H 
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D -1 (o>P>q) = i(l - i) (t f Z Z Z* dz + T 0 [ f z'z'z* dz + 
1 *• 1 J P q J 2 L J P q j 


v - 


1 I z* dzl } 
J p q j J i 


D 2 (3,P,q) = i(l + i) {t x J Z p Z* Z J dz + T 2 [ J Z p (Z*) 'z* dz 


p q j 


t- 

— f z"z*Z* dzl } 
J p q .i J J 


v - c 

D,(j,p,q) = i(l + i) {t f Z*Z Z* dz + T T f (Z*) Vz* dz + 
3 ^ 1 J P q J 2 L J ' p' q j 


Y - 


- 1 < z ? "Vi dz ] } 


D u (j,P,q) = 1(1 - i) {t f Z*Z*Z* dz + T ? [ f (Z*) '(Z*) ' Z* dz 
4 Lljpqj 2 L J p' ' q' j 


P q 2 ) 


Y - 


± J (Z*) ”Z*Z* dz 


where 


2tt 


2rr 


0 

1 


p q o 


]} 


T, = I 0 0 0. d9 f r'r 'R. rdr + [ @'®'©. d9 f R R R ^ 1 — s 2 (p) 

1 JpqJ J P q J J p q j J p q j r 2 mrr^' 

00 00 


2tt 


2tt 


X [ 0 0 ®. d9 f R R R.rdr 

J p qj J pqj 
0 


T = ©00. d9 f R R R.rdr 

2 Jpqj Jpqa 


(c-5) 


(c-6) 


(C-7) 


(C-8) 


x 
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In the equations on the prior page the notation of Eq. ( 9 ) is used; that is, 
a single index (i.e., j, p, or q) is used to identify a particular series term 
rather than the mode numbers used in Eq. ( 6 ). The index j identifies the 
equations in which a given coefficient appears which corresponds to the weight- 
ing function used in deriving that equation. For the coefficients of the linear 
terms (i.e., the C's) the index p identifies the amplitude function which the 
coefficient multiplies. For coefficients of the nonlinear terms, (i.e., the D’s) 
p identifies the factor which is not differentiated with respect to time, (i.e., 
A p or Ap ), while q identifies the differentiated factor (i.e. dAp/dt or dAp*/dt . 
Due to the complex nature of the axial eigenfunctions, the above coefficients 
are complex numbers . 

Structure of the Numerical Calculations 

A flow chart for Program C0EFFS3D is shown in Figure (C-l) . The program 
can be divided into five major sections: (l) input, ( 2 ) calculation of the 
complex linear coefficients, ( 3 ) calculation of the complex nonlinear coeffi- 
cients, (4) obtaining coefficients of the equivalent uncoupled real system, 
and ( 5 ) output . 

The inputs to the program include the various parameters describing 
the chamber geometry, the nozzle boundary condition, the modes included in the 
approximating series expansion, and various control numbers, as well as the 
roots of the Bessel functions. 

In the second section the axial acoustic eigenvalues are calculated by 
means of Subroutines EIGVAL and FCNS, and the integrals of the products of 
two axial eigenfunctions are computed by means of Subroutines AXIAL1 and UBAR. 

The integrals involving radial and tangential eigenfunctions are evaluated by 
using the orthogonality properties of these functions. The complex linear co- 
efficients are then calculated according to Eqs. (C-l) through ( C-b ) and are 
normalized by dividing by C Q (j, j). 

In the third section the integrals of products of three Bessel functions 
are calculated using Subroutines RADIAL and JBES, while similar integrals in- 
volving azimuthal eigenfunctions and axial eigenfunctions are computed using 
Subroutines AZIMTL and AXIAL2 respectively. The normalized complex nonlinear 
coefficients are obtained from Eqs. (C-5) through (C- 8 ) by dividing by C Q (j, j). 

In the fourth section the normalized complex coefficients are used to 
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obtain the coefficients for the equivalent system of real differential equations 
obtained by separating the real and imaginary parts of the complex equations . 
Since the axial eigenfunctions are non -orthogonal, the resulting system of 
equations may be coupled in the second derivative terms. Therefore, a matrix 
inversion procedure is used to obtain the coefficients of an equivalent system 
which is not coupled in the second derivatives. 

In the last section the computed values of the coefficients are either 
printed out, punched onto cards, or stored on drum (FASTRAND file) as desired. 
Input Data 

The input data consists of the chamber parameters (i.e., ratio of speci- 
fic heats, steady state Mach number, and length-to-diameter ratio), the nozzle 
admittance ratio, various control numbers, and information indicating which 
modes are included in the approximate series expansion. Regarding the latter 
information, each term in the series is identified by the integer variable J. 

The nature of each term is specified by the four integers L(j) , M (j) , N(j), 
and NS(j), and each term is given a four character name NAME(j) . In this 
manner the coefficients are identified by the integers J associated with the 

modes involved rather than the corresponding axial, azimuthal, and radial 
mode numbers . 

The followi ng comments pertain to the detailed description of the input . 
The location number refers to columns of the card. Three formats are used for 
input: "A" indicates alphanumeric characters, "I" indicates integers, and 

F indicates real numbers with a decimal point. For the "I" and "F" formats 
the values are placed in fields of five and ten locations, respectively, and 
the numbers must be placed in the rightmost locations of the allocated field. 


No . of 


Cards 

Location 


Input Item 

Comments 

1 

1-72 

A 

TITLE 

Title of Case 

1 

1-10 

F 

GAMMA 

Ratio of specific heats, y . 


11-20 

F 

UE 

Steady state Mach number at 
nozzle entrance, u . 


21-30 

F 

RLD 

ti 

Length-to-diameter ratio, 
L/D = Zg /2 


31-40 

F 

ZCOMB 

Lenth of combustion zone, 


z / 
c/z . 
' e # 
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No. of 




Cards Location 

Type 

Input Item 

Comments 

41-45 

I 

NDR0PS 

If 0: droplet momentum 
source neglected. If 1: 
droplet momentum source 
included . 

46-50 

I 

NOZZLE 

If 0: quasi-steady nozzle. r 

If 1: conventional nozzle. 

1 1-5 

I 

NJMAX 

Number of series terms (com- t 
plex) . (NJMAX a: 10) 

6-10 

I 

NONLIN 

If 0: linear terms only. 

If 1: linear and nonlinear 




terms . 

11-15 

I 

NEGL 

If 0: Nonzero coefficients 

calculated . 

If 1: Small coefficients 

neglected. 

16-20 

I 

NOUT 

If 0: printed output only. 

If 1: printed and written 




into FASTRAND file. 

If 2: FASTRAND only. 

If 3: card output only. 

If NEGL = 1: 




1 1-10 

F 

SMI 

Linear coefficients with 
absolute value less than 
SMI neglected. 

11-20 

F 

SM2 

Nonlinear coefficients with 
absolute value less than 
SM2 neglected. 

End of input for NEGL = 

1 . 



If NOZZLE = 1: 




NJMAX 1-5 

I 

J 

Integer which identifies 
series term. 

6-15 

F 

AMPL(J) 

Amplitude factor of nozzle 
admittance, A. 

16-25 

F 

PHASE (j) 

Phase of nozzle admittance, cp. 

End of Input for NOZZLE 

= 1 . 


NJMAX 1-5 

I 

J 

Integer which identifies 
series term. 

6-10 

I 

L(J) 

Axial mode number, K>. 
(0sl,(j)a£L0) 
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No. of 



Cards 

Location 


Input Item 

Comments 



11-15 

I 

M(J) 

Tangential mode number , 
(0£L(J)<;8) 



16-20 

I 

N(J) 

Radial mode number, n. 
(0sai(j)*5) 



21-25 

I 

NS(j) 

NS(j) = 1: 0. = sin(m9) 

NS(j) = 2: 0^ = cos(m9) 

* 


26-30 

A 

name(j) 

tJ 

Four character name. 


The first card gives a title (maximum 72 characters) used to identify 
the run. The second card gives the chamber parameters (i.e., y, u e , L/D, z q ) , 
determines whether the droplet momentum source is included in the analysis 
(see Appendix A) , and specifies the type of nozzle (quasi-steady or conven- 
tional) . If a quasi -steady nozzle is specified the nozzle admittance is cal- 
culated using Eqs . (l4) , and no further information concerning the nozzle is 
required. The contol numbers are given on the third card. Due to computer 
storage limitations the series expansion is limited to ten terms, thus NJMAX ^ 
10. The control number NEGL gives the option to neglect all coefficients 
with absolute value smaller than a given number, thus allowing a considerable 
saving in computation time when the equations are numerically integrated by 
Program LCYC3D. It has been found that neglecting coefficients with absolute 
value smaller than 0.1 (i.e., SMI = SM2 = O.l) reduces the computation time 
by half and has a negligible effect on the resulting solutions . For conven- 
tional nozzles a series of NJMAX cards is read which gives the nozzle ad- 
mittance (amplitude and phase) for each term in the series. This is followed 
by another series of NJMAX cards giving the mode numbers for each series term. 

The proper input for program C0EFFS3D will be illustrated with the 
following example. Suppose the velocity potential § is expressed in terms of 
the first tangential (IT), the second tangential (2T), and the first radial 
(1R) modes. It is also desired to investigate instability of the spinning 
type, therefore both sin(me) and cos(m9) terms are included in the series. 
However, for the 1R mode (m=0) there is no corresponding sin(m9) term, there- 
fore the resulting series will contain five terms. A nozzle admittance of 
A = 0.02 and cp = 45 will be assumed for each term in the series, and coeffic- 
ients smaller than 0.1 as well as the droplet momentum source will be neglected. 
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The output data will be punched on cards. A sample input for this case is 
given in Table (C-l) below. 


Table C-l. Sample Input. 


Q 

IQ 

1! 

3G 

■u 

1 

IQ 

IQ 

■ 

mu 

IS 

IS 

mu 

IB 

mu 

IQ 

HI 

Ifl 

HI 

IB 

mu 

Q 

mu 

0 

mu 

II 

HI 

II 

mu 

I 

mu 

II 

tj 

■ 

u 

1 

1 

u 

1 

KJ 

1 

mu 

1 

til 

1 

u 

1 

1 

til 

1 

til 

1 

u 

1 

0 

1 

a 

1 

El 

1 

EJ 

■ 

El 

■ 

EJ 

1 

EJ 

■ 

0 

1 

0 

1 

ra 

■ 

□ 

■ 

KU 

1 

Cj 

1 

0 

1 

E3 

1 

ill 

B 

n 


1 

ii 

III 

II 

■ 

■ 

Q 

II 

10 

II 

II 

1 

1 

I 

II 

1 

0 

1 

0 

1 

i 

1 

1 

1 

1 

1 

0 

1 

0 

1 

1 

l 

1 

1 

1 

■ 

0 

1 

0 

1 




0 

I 


I 

| 

91 

B 

m 


» V » 

w to 

8 

B 

8 

K9 

n 

Bl 

Bl 

Bl 

8 



i 

III 

IS 

II 

II 


II 

IQ 

II 

II 

II 

II 

IQ 

II 

II 

1 

1 


1 

i 

1 

1 

1 

1 

1 

1 

1 

1 

1 

l 

l 

1 

1 

| 

1 

1 

| 

I 

i 

■ 

■ 

■ 

■ 

III 

■ 

gg 

B 

m 


n 

Bin 

8 

IS 

D 

B 

K9 

19 

K9 


KS 


1 

■ 

II 

II 



0 

1 

a 

■ 

I 

I 

I 

1 

1 

1 

0 

1 

Q 

r 




1 



" 

" 


r 

" 

" 

" 

1 

r 

r 

r 

r 

1 

1 

1 

1 

1 

1 

| 


I 

I 

n 

E5 

KB 

IB 

n 

n 


» V It Jf 


8 

19 


8 

8 

m 

8 

m 

El 


u 

II 

IQ 

1 

II 

1 

I 

1 

1 

10 

1 

10 

10 

II 

I 

II 

I 

i 

ID 

0 

I 

0 

■ 

■ 

1 

1 

■ 

■ 

■ 

I 

1 

1 

II 

i 

■ 

■ 

■ 

1 

I 

■ 

■ 

1 

ill 

1 

H 

m 

D 

KB 


a 

m 

8 


S3 


1 

i 

II 

IQ 

1 

1 

1 

1 

1 

■ 

0 

1 

0 

0 

1 

1 

1 

■ 

1 

i 

B 

5 

1 

0 

1 

1 

1 

1 

1 

i 

1 

1 

1 

■ 

1 

i 

■ 

i 

■ 

■ 

1 

■ 

■ 

1 

II 

l 

I 

II 

BM 


8 


B 

B 


III! 

is 

1 

1 

1 

■ 

1 

1 

0 

1 

0 

0 

1 

1 

1 

1 

1 

i 

9 

9 

1 

0 

1 

1 

1 

1 

1 

i 

1 

1 

1 

1 

1 

i 

■ 

i 

■ 

1 

1 

■ 

B 

1 

l 


1 

1 

II 

D 

n 

KB 


a 

B 

B 

a 

B 

B 





ID 

1 

1 

1 

■ 

1 

1 

0 

1 

0 

0 

1 

1 

1 

1 

1 

i 

B 

9 

1 

0 

1 

1 

1 

1 

1 

i 

1 

1 

1 

1 

1 

i 

i 

i 

■ 

1 

1 

■ 

■ 

1 

■ 

i 

1 

1 

ii 

n 

n 

Bfl 


B 

B 

B 



1 

II 

la 

1 

1 

1 

■ 

1 

■ 

0 

B 

0 

0 

1 

1 

1 

■ 

1 

i 

S 

9 

1 

& 




1 

1 

i 

1 

1 

1 

1 

■ 

i 

■ 

i 

■ 

1 

1 

■ 

■ 

1 


1 

1 



n 



8 

8 

PI 

B! 



III 

IQ 

1 

1 



0 

1 

1 

1 

1 

0 

1 

1 

1 

1 

Q 

i 

1 

1 

1 

Q 

1 

□ 

0 

Q 

9 

r 



" 

' 


" 

n 


i 

“ 

" 

n 

n 

n 

□ 

n 

III! 

n 

n 

KB 

Bl 

n 

n 

n 

a 

Bl 


4* 

if * *9 JO W 
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1 

1! 

■ 

II 

i 
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1 

1 

1 

0 

1 

1 

1 

1 

Q 

i 

1 

1 

1 

0 

1 

0 

0 

0 

0 











" 

“ 

n 


- 



“ 

n 

T1 

n 

n 

KB 


B 

B 

B 

B 

B 
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i 

II 

II 

IS 


1 

1 


0 

1 

1 

1 

1 

0 

1 

1 

■ 

1 

0 

i 

1 

1 

1 

0 

1 

Q 

0 

0 

9 

1 

1 

1 

1 

1 

1 

i 

■ 

i 

■ 

1 

1 

1 

■ 

1 

1 

1 

1 

1 

■1 

n 

n 

KB 



i 

II 

ID 

1 

1 

1 

1 

0 

1 

1 

1 

1 

0 

1 

1 

1 

1 

0 

i 

1 

1 

1 

0 

1 

0 

Q 

0 

0 

1 

1 

1 

1 

1 

1 

i 

■ 

i 

■ 

1 

1 

1 

1 
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1 


n 
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Bl 
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II 
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■ 

II 


01 


II 



0 


1 



0 


1 

1 

1 

0 

1 

9 

Q 

0 

0 

1 

1 

■ 

■ 

1 

1 

i 

i 

i 

i 

: 

: 

: 

□ 


□ 

□ 

□ 

□ 

u 


After the last card in the sequence described above is read, the program 
is executed and control returns to the input section. Thus, several cases can 
be executed on the same run. If no further cards are given the run is termi- 
nated. 

In addition to the above card input, roots of the Bessel functions S 
, . , . mn 

which give zero slope at r - 1 and the associated values J (S ) are needed for 

m mn 

these calculations. These values were taken from Ref. (l6) for m = 0, 1,...8 
and n - 1,2... 5; they are automatically put into the program by means of a DATA 
statement, which is an integral part of the program. 
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Complex Linear Coefficients 


For NDROPS = 0 the complex linear coefficients are computed from Eqs. 
(C-l) through (C-4) and are stored in the complex array CC(KC,NJ,NP) . For 
NDROPS = 1 the coefficients C 2 (j,p) are computed from Eq. (A-8). 

In order to calculate these coefficients the following information is 
needed: (l) the axial acoustic eigenvalues, b^, (2) the steady state Mach 

number distribution, u(z), (3) the orthogonality properties of the trans- 
verse eigenfunctions, and (4) the integrals of products of two axial eigen- 
functions. The calculation of these quantities is described below. 

Axial Acoustic Eigenvalues . The axial acoustic eigenvalues are deter- 
mined by numerically solving the transcendental equation given by Eq. (8) . 
This is done by first substituting b^ = + iT)^ and Y = + iY i 

into Eq. (8) and separating real and imaginary parts. This yields a pair 
of simultaneous equations of the form: 


f( e j7]) = 0 
g(e,T|) = 0 


(c-9) 


where 


f(e.ll) - (e 2 - T| 2 )f( 6)T|) - 4«T|H(«,T|) 

+ y 2 {[ (Y r - Y i><i + ' 2 - - '•Vi'll <K«,1l) 

+ “CwC + s 2 - n 2 ) + -A} ( C . 10) 
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g(e?T]) - (e -7) )H( e ,T]) + eT)F( € >T)) 


+ ^{[W 3 ™ + s' - D 2 ) + < Y r - 

- [ {y2 r - Y i>< S L + = 2 ' /) + ( C - Y 

and 

F( e ,Tl) = sin 2 ( ez e )cosh 2 (i]z e ) - cos 2 (ez e )sinh 2 (7)z e ) 

G(e,T|) = cos 2 (ez e )cosh 2 (T)z e ) - sin 2 ( e Ze )sinh 2 (T]z e ) (C-l 

H(e>T)) = sin(ez )cos(ez )sinh( 7 ]z )cosh(Tlz ) 

c c 0 G 


In the above equations the subscripts on e and 7] have been omitted. 

Equations (C-9) are solved by Subroutine EIGVAL using Newton’s Method 

17 

for two unknowns. In this method successive approximations to the roots 
are generated by the recursion formulas: 


e i+l = e i - 

ri g y - fin 

L J(f,g) Ji 

= \ - 

r gf e " fg p1 

~(f,g) "Ji 


where the Jacobian J(f,g) is given by: 



and the subscripts indicate partial differentiation with respect to g and j). 
The quantities f, g, f^, f^, g g , g^ are calculated by the Subroutine FCNS. 
The iteration is started by assuming the following values for e and T): 


e o = ®m + a cos ^) 

T] o = a sin(p) 

where for t = 0: e = 0 

m 

a = 10A/ z g 

3 =9/2+45 (degrees) 

and for £ / 0: ^ = tn/z e 

a = A/z g 

3 = cp + 90 (degrees) 


(C-15) 


(C-l6) 


The iteration is terminated when the errors Ae and ^T| are smaller than lO - ^. 
If the iteration fails to converge after 40 iterations or the Jacobian 
vanishes a warning message is printed. FORTRAN listings of Subroutines 
EIGVAL and FCNS are given at the end of this appendix. 

Steady State Mach Number Distribution . The steady state Mach number 
distribution is calculated by means of Subroutine UBAR which must be sup- 
plied by the user . This distribution must be of the form shown in Fig. 

(C-2) where the Mach number varies from zero at the injector face (z = 0) 

to its maximum value at the end of the combustion zone (z = z ) and remains 

c 

constant until the nozzle entrance (z = z ) is reached. Thus the Mach 
number is given by 


u(z) = U(z)u (0 £ z ^ z ) 

e c 

(C-17) 
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u(z) = u e (z c <; z «; z e ) 

where U(0) = 0 and U(z c ) = 1. Although the function U(z) may be arbitrary, 
the results presented in this report were obtained using a linear Mach num- 
ber distribution in the combustion zone (i.e., uniformly distributed com- 
bustion). Thus the function U(z) in the listing of UBAR provided herein is 
given by: 


U(z) = z/z c (C-l8) 

In addition to the Mach number distribution (NOPT = l), the first (NOPT = 2) 
and second (NOPT = 3) derivatives are also calculated. 



Axial Coordinate , z 


Figure C-2. Steady-State Mach Number Distribution. 
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Orthogonality of Transverse Eigenfunctions . The tangential eigen- 
functions have the following orthogonality properties: 


r 2jT 

i 0 s 


r . 2 TT 


i. 2 TT 


P 

cos(m 0)cos(m.0)de = 
J o P 3 

p2rr 

sin(m e)cos(m e)de 


= 0 

m 

m. 



P 

r 1 


= TT 

m 

= m. 



P 

J 


= 2rr 

m 

= m. 

= 0 


P 

3 


= 0 

for 

all 

m < 


(c-19) 


For the special case of m - m - 0 the integral involving sines vanishes. 

The orthogonality property of the radial eigenfunctions is given by: 

[r R.rdr = 0 n ^ n. (m = m.) 

J 0 p 3 p r j p r 

o s 2 p (c - 20) 

J/pY^ ’ [ J m< W] "p - <-p - ) 

the tangential integrals vanish when m^ ^ m. it is not necessary to 

calculate the radial integrals for m^ ^ nn . These orthogonality properties 

are used to calculate the integrals, P 11 0 I 0.de and p R R.rdr, which appear 

• ip , . , > Jupj J o p j 

in Eqs . (C-1J through (C-4). For a series containing pure transverse modes 

only (l = 0), it is easily seen that all of the linear coefficients vanish 

except those corresponding to p = j, yielding a system of equations which 

are not coupled in the linear terms. 

Axial Integrals . The integrals of products of two axial eigenfunctions 
are calculated by Subroutine AXIAL1. According to the value of the input 
parameter NOPT these integrals are calculated as follows: 
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NOPT =1: J 

, Z e , sinhfi(b + b.)z 

1 „ „* if . L P . , j/_ e 

Li L> UZ “ o'! 

0 P 3 H 

i(b + b*) 
P 3 


sinh 

L > 

fop - b >a] •. 



Kb - b!) > 

p j 

NOPT = 2: 

z 

I* e " * 2 

Z Z.dz = -b 
l 0 P 0 P . 

z 

f Z Z.dz 
’o P 

NOPT = 3: 

% 

Z 

r 6 du * 

^ Z Z.dz 
) 0 dz P J 

(evaluated numerically) 

NOPT = 4 t 

« 

z 

r e -/ \ ' * 

u(z)Z Z .dz 
J v P 3 

(evaluated numerically) 


o 


(C-21) 




(C-22) 


The last two integrals, which involve the mean flow Mach number, are eval- 
uated by means of Simpson's Rule. A FORTRAN listing of AXIAL1 is provided 
at the end of this appendix. 


Complex Nonlinear Coefficients. 

The complex nonlinear coefficients are calculated from Eqs. (C-5) 
through (C-8) and are stored in the complex arrays, CDl(NJ,NP,NQ) , 
CD2(NJ,NP,NQ) , CD3(NJ,NP,NQ), and CD4(NJ,NP,NQ) . 

In order to calculate these coefficients, the various integrals of 
axial, azimuthal, and radial eigenfunctions must be evaluated. Since many 
of the azimuthal integrals are zero they are evaluated first, and the re- 
maining integrals are computed only if the corresponding azimuthal integral 
is nonzero. The subroutines used to calculate these integrals are described 
in the following paragraphs. 

Azimuthal Integrals . The azimuthal integrals are calculated by Sub- 
routine AZIMTL according to the value of NOPT as follows: 


NOPT = 1 : 


r a,T 

J Wi d| 
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NOPT = 2 : 



/ / 

© © © .d.0 
PTJ 


These integrals are easily evaluated analytically; for most values of p, q, 
and 3 they are zero. The nonzero integrals are readily expressed in terms 
of the following integrals: 


r 


i2tt 


cos(m e)cos(m 9)cos(m.e)d0 = Tr/2 for m. = m + m , 

p q 3 3 p i 


m = 


m. + m , or 

3 q 


m = m. + m 

q j p 


(C- 23 ) 


r 


,2tt 


cos(m^e)sin(m^9)sin(m^.9)de 


tt/ 2 for m = m + m. or 

' q p 3 


m. = m + m 
3 p q 


(C- 24 ) 


/ 


,2n 


cos(m 0)sin(m 9)sin(m.9)d9 = -n /2 for m = m + m. 
P q J P q 3 


(C- 25 ) 


where m , m^, and nu are nonzero. If any one of the tangential mode numbers 
is zero (corresponding to a radial mode) the following values are obtained: 

i.2tt 



cos(m 0)cos(m 0)cos(m.0)d0 = 2tt 

j 0 p q d 

m 

P 

II 

li 

cJ 3 

= 0 

* 

= TT 

m 

P 

= o, 

11 

s 04 

V 



m 

q 

II 

O 

11 

V 



m. 

J 

= o, 

II 

m 

q 


(c-26) 
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TT 


in. 

J 


(C-27) 




m - 


Subroutine AZIMTL consists of two sections. In the first section the 

aziinuthal integral is expressed as the product of a constant factor and one 

of the basic forms given in Eqs. (C-23) and (C-24). The second section is 

essentially a series of logical tests to determine if the mode numbers, m , 

P 

m , and m satisfy any of the conditions for Eqs. (C-23) through (C-27). If 
u j 

any of these conditions is satisfied the appropriate value is mul tiplied by 
the corresponding factor determined in the first section and the product is 
assigned to the output variable (i.e., RESULT), otherwise the value zero is 
assigned. 

Radial Integrals . Subroutine RADIAL calculates the radial integrals 
which appear in Eqs. (C-5) through (C-8) according to NOPT as follows: 


NOPT = 1 : 
NOPT = 2 : 

NOPT = 3 : 


£ 

t 


R R R.rdr 
P ? J 


RRRidr 

o p 1 d- 


f^R Vr . 
J„p<n 


rdr 


where the R's are the Bessel functions, JjS^r). These integrals are com- 
puted numerically using Sinqpson's Rule with 100 subdivisions. In calculating 
the integrands the derivatives of the Bessel functions are given by: 

J m ( S mn r ) = i [ J m-l( S mn r ) ’ t Vl( S mn r )] for m = 1,2,3,... 

(C-28) 

J o( S mn r ) “ _J l( S mn r ) 
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The integrand of the second integral (NOPT = 2) is indeterminate at the lower 
limit of integration. However a limit exists, denoted by L, which vanishes 
with the following exceptions: 


L = S (p)/2 
mn ' 

for 

m = 1, 
P 

II 

II 

o 


L = 

for 

II 

H 

o 

II 

gp 

II 

(c-29) 

1 = V<3>/2 

for 

«\ 

i — I 

ii 

eP 

II 

II 

o 



% 
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All of the calculation^ in Subroutine RADIAL are carried out in double 
precision arithmetic. The results are given as a single precision number. 

Subroutine JBES computes the double precision Bessel functions which are 
needed for the above calculations. A description of this subroutine and a 
program listing are given in Chapter 23 of Ref. (l8). 

A xial Integrals . The integrals of the products of three axial eigen- 
functions (see Eqs . (C-5) through (C-8)) are computed by Subroutine AXIAL2 
according to the input parameters NOPT and NCONJ. The three basic forms 
are specified by NOPT as follows: 


NOPT = 1 ; 


NOPT = 2 : 


NOPT = 3 : 


f 


Z Z Z*dz 
py 


J 


o W* 


i» e " * 

z Z Z.dz 

J 0 p ^ J 


When NCONJ = 1 these basic forms are calculated; these are the forms appearing 
in the expression for D^j, p, q) (see Eq. (C-5)). For NCONJ = 2 the second 
function in the integrand is replaced by its complex conjugate to obtain the 
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integrals appearing in the expression for D 2 (j,p,q). The integrals appearing 
in the expressions for D 3 (j,p,q) and are obtained by setting 

NCONJ = 3 and NCONJ = 4 respectively. 

The basic forms are calculated from the following analytical formulas: 


p Z e * , r sinh|i(b + b + b*)z 1 

r z z z dz = i { L_E 9 — -J—SJ 

«L P q 3 U Wv ^ a. ^ 


i(b + b + b .) 

p q y 

sinh|”i(b + b - b*)z 1 

Ll.p a . y eJ 

i(b + b - b*) 

P q 3 

sinh[i(b - b + b*)z 1 

L E a 3 ej 

i(b - b + b*) 
sinh|i(b - b - b*)z 1 

k E a 3 eJ 


J 


- ' ' * 

Z Z Z.dz 
P q J 



i(b - b - b*) 

p q 3 

h h 

. sinh i(b + b 

/ L P q 

p 

q ^ i(b + b + 

p q 

sinh 

i(b + b - b.)z ] 
L P q 3 ej 

i(b + b - b*) 
P q 3 

sinh 

i(b - b + b*)z 1 
L P q 3 el 

l(b p - b q + b l> 

sinh 

i(b - b - b .)z 

L e q j e - 


} 






i(b - b - b.) 

p q 3 


(C-30) 


(C-31) * 


I 


e " * 

Z Z Z.dz 

. p q 3 


2 f * 

- b Z Z Z.dz 

p J Q p q 3 


(C-32) 
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The remaining forms are obtained from Eqs. (C-30) through (c-32) by replacing 

the appropriate eigenvalues with their complex conjugates; thus, for NOPT - 2 

b is replaced by b* for NOPT = 3 b is replaced with b\ and both b and 

4. P P P 

b^ are replaced by their conjugates for NOPT =4. 

FORTRAN listings for Subroutines AZIMTL , RADIAL, and AXIAL2 are given at 
the end of this appendix. 

Coefficients for Equivalent Real System . 

Equations (12) are a system of complex differential equations to be 
solved for the unknown complex amplitude functions, A (t). In order to solve 
these equations numerically they must first be separated into their real and 
imaginary parts. This is done by assuming that A p (t) = F (t) + iG (t), sub- 
stituting into Eqs . (12), and separating real and imaginary parts to obtain 
an equivalent system of real differential equations that describe the behavior 
of the F p s and G^'s. Since these equations contain twice as many unknown 
functions (i.e., F p (t) and G^(t)) as Eqs. (12), it is convenient to re-iridex 
the unknown functions and their coefficients as follows: 


F p (t) *= B 2p-1^> 


G p (t) = V t} 


(c-33) 


Thus the B's with odd indices correspond to the real parts, F (t), and the 
B's with even indices correspond to the imaginary parts, G (t). The corre 
sponding set of differential equations is given by: 


2N 


Z { c o(j>p) 

P-1 



+ C 1 (j,p)B p (t) + [C 2 (j,p) - 



+ 
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The real coefficients in Eqs. (C-34) (i.e., C^C^C^C^, and D ) are related 
to the complex coefficients in Eqs. (12) (i.e., D^,...D^) as 

follows: 

c'(2j-l, 2p-l) = Re [C k (j,p)] 

c£(2j-l, 2p) = -Im [c k (j,p)J 

(C-35) 

C k ( 2 j, 2 p-l) =lm[c k (j,p)] 

C fc ( 2j , 2p) = Re [c k (j,p)J 

for k = 0,1,2, 3, j = 1,2, ...N, p = 1,2, ...N and: 

D (2j-l,2p-l,2q-l) = Re ^(^Pjq) + D 2 (j,p,q) + D 3 (j,p,q) + ,P,q) J 

D (2j-l,2p-l,2q) - Im ^-D 1 (j,p,q) + Dg(;J,p,q) - D 3 (j,p,q) + D^(j,p,q)J 

D (2j-l,2p,2q-l) = Im ^-D 1 (j,p,q) - D 2 (j,p,q) + D 3 (j,p,q) + D 4 (j,p,q)J 

D (2j-l,2p,2q) = Re [-D-^j^q) + D 2 (j,p,q) + D 3 (j,p,q) - j ,P ,q)] 

(C-36) 


6o 


D ( 2j , 2p-l,2q-l) = 
D (2j,2p-l,2q) = 
D (2j , 2p,2q-l) 

D (2j,2p,2q) 


Iffi [V j ’ P ^) + ° 2 ( + D 3 (j 5 P>q) + D 4 (j,p,q)J 
Re [D^j.p.q) - D 2 (j,p,q) + D 3 (j,p,q) - D 4 (j,p,q)J 
Re ^D 1 Cj,p,q.) + D 2 (j,p,q) - D 3 (j,p,q) - D 4 (j,p,q)J 

Im [-V^’P’q) + + d 3 (j>p><i) - D 4 (j,p,q)] 


for j 1,2, .. .N, p - 1,2,...N, q = 1,2,...N. The linear coefficients 
are stored in the arrays Cl(NJ,NP) for k = 0 and C(KC,NJ,NP) for k = 1,2,3. 
The nonlinear coefficients are stored in the array D(NJ,NP,NQ) . 

In general Eqs. (C-34) are coupled in the second derivatives; that is, 
they are of the form: 


2N / d 2 B 

I { C 0<J.P) —£ } = g 3 (B 1 ,B 2> ... B 2K ) (C-37) 

P = 1 Z 

/ 

where there are two or more C Q terms in each equation. This coupling results 
from the non-orthogonality of the axial eigenfunctions. In order to numeri- 
cally integrate Eqs. (C-34), they must be decoupled by transforming to the 
form: 


d 2 B 


dt 




T 1’ 2 


(c-38) 


in which only one second derivative appears in each equation. Using Eq. 
(C-38), it is seen that Eq. (C-37) can be expressed as 


c 0 f = « (c-39) 

/ 

where Cq is the 2N X 2N matrix of coefficients of the coupled system, f is 


6l 


the column matrix corresponding to the right-hand- side of the decoupled 
system, and g is the column matrix corresponding to the right -hand -side of 
the coupled system. To decouple Eqs. (C-37), therefore, Eq. (C-39) is solved 
for f, thus: 

f = C^g (C-40) 

_1 > 

where C Q is the inverse of the matrix C Q . Performing these operations and 
equating the coefficients of like terms in f and C^g gives the following 
relations: 


2N 

c i(j> p ) =£ C~ 1 ( t j ,k)c!(k,p) i = 1,2,3 
k=l 

2N 

D(j,P,q) = ^CQ 1 (j,k)D (k,p,q) 
k=l 


(C-41) 


/-s»/ 

where C. and D are the corresponding coefficients of the decoupled system. 

1 _1 

The matrix inverse, C Q , is computed by the subroutine GJR, which is a 
standard Univac 1108 library program, and is stored in the array Cl(NJ,NP). 

A listing of GJR and instructions for its use are given in Ref. (19)* 

rsj r>*» ' 

The calculation of C^(j,p) and D(j,p,q), which are the coefficients for 
the equivalent set of real, decoupled equations, is the final step in the 
computations performed by C0EFFS3D. The coefficients are stored in the arrays 
C(KC,NJ,NP) and D(NJ,NP,NQ), replacing those computed from Eqs. (C-35) and 
(C- 36 ) . The output of these coefficients is described below. 

Output 

According to the value of the control number NOUT, the coefficients 
calculated by Program C0EFFS3D are printed, punched onto cards, or stored on 
drum (FASTRAND). These three output modes will now be discussed indivi- 


dually. 

Printed Output . Since the printed output cannot be used as input to 
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Program LCYC3D, the option "printed output only" (NOUT = 0) is only used for 
checkout purposes. Printed output can also be obtained in conjunction with 
the drum storage mode (NOUT = l) . Since the printed output format can only 
accommodate five series terms (complex), it should only be used for NJMAX £ 5 . 

The first page of printed output gives a restatement of the input para- 
meters. This page is headed by the title of the case (TITLE) which is fol- 
lowed by the ratio of specific heats (GAMMA), the steady state Mach number 
at the nozzle entrance (UE), the length-to-diameter ratio (L/D), and the 
length of the combustion zone as a fraction of the chamber length (ZCOMB) . 
After statements concerning the presence or absence of the liquid droplet 
momentum source and the type of nozzle considered, a restatement of the input 
parameters J, L(j), M(j), N(j), NS(j), and NAME(j) which describe the terms 
in the series expansion of $ is given. This tabulation also includes addi- 
tional parameters needed by Program LCYC3D: S , the dimensionless frequen- 

cy of the mode (SMN): J (S ), the associated value of the Bessel function 
(JM(SMN)); the real part (EPS) and the imaginary part (ETA) of the axial 
acoustic eigenvalue; and the real part (YR) and imaginary part (Yl) of the 
nozzle admittance. 

The next three pages give the decoupled linear coefficients, CL(j,p), 
Cp(jjp)} and C^jjp). These coefficients are presented in the matrix format 
with the rows corresponding to the index j and the columns corresponding to 
the index p. The remaining pages give the decoupled nonlinear coefficients 
U(j,p,q) for each value of j . Here the rows correspond to the index p and 
the columns correspond to the index q. 

A sample printed output for the five term series used in the sample 
input is given in Tables (C-2) through (C-4). 

Drum Storage . When available drum storage, such as the FASTRAND system 
used with the Univac 1108, is the most convenient means of storing the out- 
put of Program C0EFFS3D. In the absence of such a system, the program can 
be easily modified to store the coefficients on magnetic tape. In either 
case magnetic tape can be used as a back-up file or for permanent storage 
of the data. The control statements needed to execute these procedures 
depend upon the computer facilities being used and cannot be described in 
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Table C-2. Sample Printed Output, Page 



6k 


Table C-3. Sample Printed Output, Page 
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Table C-4. Sample Printed Output, Pa 
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QOOOOO* 000000* 000000* 000000* oonooo* 000000* 000000* 000000* 968010*- T 1©296* 


this manual. 


Card Output . When a drum or magnetic tape storage is not available, 
punched card output can be used (NOUT = 3). This method becomes unwieldy, 
however, when a large number of coefficients is involved since only one co- 
efficient can be punched on a card. The format for both drum and card out- 
put is the same and is given below: 


Number 





of Cards 

Location 

Type 

Output Item 

Comments 

1 

1-10 

F 

GAMMA 

Same as for input. 


n-20 

F 

UE 

Same as for input. 


21-30 

F 

ZE 

Dimensionless chamber length, 
(2L/D). 


31-40 

F 

ZCOMB 

Same as for input. 


4l-45 

I 

NDROPS 

Same as for input. 


46-50 

I 

NJMAX 

Number of unknown functions, 
B p (t) (see Eq. (C-34)). 

njmax/2 

1-5 

I 

J 

Same as input. 


6-10 

I 

L(J) 

IT 


11-15 

I 

M(J) 

II 


16-20 

I 

N(J) 

II 


21-25 

I 

NS( J) 

IT 


26-35 

F 

S(J) 

Root of Bessel function. S 

5 mn 


36-45 

F 

SJ(J) 

Associated value of Bessel 

function, J (S ). 

nr mn 


46-50 

A 

NAME(J) 

Same as input. 

NJMAX/2 

1-5 

I 

J 

Same as input. 


6-15 

F 

YR 

Real part of nozzle admit- 
tance, Y . 

r 


16-25 

F 

YI 

Imaginary part of nozzle 

admittance, Y. . 

1 


26-35 

F 

EPS 

Real part of axial eigen- 
value , e . 


36-45 

F 

ETA 

Imaginary part of axial 
eigenvalue, 7 ]. 
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Number 


of Cards 

Location 

T£E£ 

Output Item 

Comments 

1 

1-5 

I 

KMAX(l) 

Number of nonzero linear 
coefficients of type O^Cjjp) 

KMAX(l) 

1-5 

I 

NJ 

Index, j. 


6-10 

I 

NP 

Index, p. 


11-25 

F 

C(l,NJ,NP) 

Linear coefficient, ,p) . 

1 

1-5 

I 

KMAX(2) 

Number of nonzero linear 
coefficients of type C 2 («j,p) 

KMAX(2) 

1-5 

I 

NJ 

Index, j . 


6-10 

I 

NP 

Index, p. 


11-25 

F 

C(2,NJ,NP) 

Linear coefficient, C 2 (j ,p) . 

1 

1-5 

I 

KMAX(3) 

Number of nonzero linear 
coefficients of type C^(j,p) 

KMAX(3) 

1-5 

I 

NJ 

Index, j . 


6-10 

I 

NP 

Index, p. 


11-25 

F 

C(3,NJ,NP) 

/N/ 

Linear coefficient, C^( j , p) . 

1 

1-5 

I 

kmax(4) 

Number of nonzero nonlinear 
coefficients . 

KMAX(4) 

1-5 

I 

NJ 

Index, j . 


6-10 

I 

NP 

Index, p. 


11-15 

I 

NQ 

Index, q. 


16-30 

F 

D(NJ,NP,NQ) 

Nonlinear coefficient, 
D( • 


The first card of output gives the chamber parameters y, u^, L/D, and 
z c / z e ’ the dr °P le ' t momentum source control number, NDROPS; and the n umb er 
of unknown real functions (i.e., B (t)), NJMAX. This is followed by 
NJMAX/2 cards (the number of unknown complex functions, A (t)) describing 
the terms included in the series expansion of $. The next NJMAX/2 cards 
gives the complex nozzle admittance (Y^ and Y^) and the corresponding com- 
plex axial eigenvalue (e and j|) for each complex series term. The linear 
coefficients are given in three sets of cards. The first card in the set 
gives the number of coefficients of the given type, while the remaining 
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cards give the indices j and p and the coefficient Ch(j,p). The next card 
gives the number of nonlinear coefficients and is followed by cards giving 
the indices j, p, q and the corresponding coefficient D(j,p,q). Both linear 
and nonlinear coefficients are given in a field of 15 spaces with six deci- 
mal places. For NEGL = 0 only the nonzero coefficients (absolute value 
greater than 10 are given, while for NEGL = 1 only linear coefficients 
with absolute value greater than SMI and nonlinear coefficients with abso- 
lute value greater than SM2 are given. 

A sample card output produced by the sample input of Table (C-l) is 
given in Table (C-5) below. 


Table C-5* Sample Card Output. 


> i * >6? » t ion u n i< n it i7 » n 79 ji » n u » w v » » v i? inontrii n «o o ** *} 
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FORTRAN Listing . 


********************* PH0GRAM C0EFFS3D *************************** 


THIS PROGRAM COMPUTES THE COEFFICIENTS WHICH APPEAR 
IN THE DIFFERENTIAL EQUATIONS WHICH GOVERN THE MODE-AMPLI TUDE 
FUNCTIONS. THESE COEFFICIENTS ARE STORED ON DRUM OR 
PUNCHED ONTO CARDS FOR INPUT INTO PROGRAM LCYC3D* 

THE FOLLOWING INPUTS ARE REQUIRED: 

THE TITLE OF THE CASE. 

GAMMA IS THE SPECIFIC HEAT RATIO. 

UE IS THE STEADY STATE MACH NUMBER AT THE NOZZLE ENTRANCE. 

RLD IS THE LENGTH- TO -DIAMETER RATIO. 

ZCOMB IS THE LENGTH OF THE REGION OF UNIFORMLY DISTRIBUTED 
COMBUSTION* EXPRESSED AS A FRACTION OF THE CHAMBER LENGTH. 
NDROPS DETERMINES THE PRESENCE OF DROPLET MOMENTUM SOURCES: 
NDROPS = 0 DROPLET MOMENTUM SOURCE NEGLECTED. 

NDROPS = 1 DROPLET MOMENTUM SOURCE INCLUDED. 

NOZZLE SPECIFIES THE TYPE OF NOZZLE USED: 

NOZZLE = 0 QUASI -STEADY. 

NOZZLE = 1 CONVENTIONAL NOZZLE. 

FOR CONVENTIONAL NOZZLE: 

AM PL IS THE NOZZLE AMPLITUDE RATIO. 

PHASE IS THE NOZZLE PHASE SHIFT. 

NJMAX IS THE NUMBER OF MODE-AMFLI TUDE FUNCTIONS IN THE ASSUMED 
SERIES SOLUTION. NJMAX MUST NOT EXCEED 10. 

THE COEFFICIENTS COMPUTED ARE DETERMINED BY NONLIN AS FOLLOWS: 
NONLIN = 0 LINEAR COEFFICIENTS ONLY. 

NONLIN = 1 BOTH LINEAR AND NONLINEAR COEFFICIENTS. 
COEFFICIENTS TO BE NEGLECTED ARE DETERMINED BY NEGL 
AS FOLLOWS: 

NEGL = 0 TERMS SMALLER THAN 0.00001 ARE NEGLECTED. 

NEGL = 1 LINEAR TERMS SMALLER THAN SMI AND NONLINEAR 
TERMS SMALLER THAN SM2 ARE NEGLECTED. 

THE OUTPUT IS DETERMINED BY NOUT AS FOLLOWS: 

NOUT = 0 PRINTED OUTPUT ONLY. 

NOUT = 1 PRINTED AND STORED ON DRUM (FASTRAND FILE). 

NOUT = 2 FASTRAND FILE ONLY. 

NOUT = 3 CARD OUTPUT ONLY. 

EACH MODE-AMFLI TUDE IS ASSIGNED AN INTEGER J. 

THE MODE IS SPECIFIED BY THE INDICES L(J)* MCJ)* AND N(J). 

LCJ) IS THE AXIAL MODE NUMBER AND MUST NOT EXCEED 10. 

M(J) IS THE AZIMUTHAL MODE NUMBER AND MUST NOT EXCEED 8. 

N(J) IS THE RADIAL MODE NUMBER AND MUST NOT EXCEED 5. 

THE INTEGER NSCJ) IS ASSIGNED AS FOLLOWS: 

NS « 1 A-FUNCTI ON SIN(M* THETA) * COSH<I*B*Z) 

NS e 2 B-FUNCTION COSCM* THETA) * COSH(I*B*Z) 

NAME(J) ISA FOUR-CHARACTER NAME. 
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c 

c 


c 

c 

c 


c 

c 


c 

c 


**************************************************»**,,********»* 
DIMENSION 

1 

2 

3 

4 

COMPLEX 

1 
2 

3 

4 

5 

COMMON 

DATA INPUT. 

PI « 3.1415927 
SMI « O.OOOOl 
SM2 ■ 0.00001 
Cl * (0*0. 1.0) 

INPUT BOOTS AND VALUES OF BESSEL FUNCTIONS. 

DATA ( ( RJROO T< I . J ) . J * 1.5). I = 1.9)/ 

1 3.83171. 7*01559. 10.17347. 13.32369. 16.47063. 

2 1.84118. 5.33144. 8.53632. 11.70600. 14.86359. 

3 3*05424. 6*70613. 9.96947. 13*17037. 16.34752. 

4 4.20119 . 8.0 1 524. 11.34592. 14 . 58 585. 17.78 8 7 5. 

5 5.31755. 9.28240. 12.68191. 15.96411. 19.19603. 

6 6*41562. 10*51986. 13*98719. 17*31284. 20*57551. 

7 7.50127. 11.73494. 15.26818. 18.63744. 21.93172. 

8 8.57784. 12.93239. 16.52937. 19.94185. 23.26805. 

9 9*64742. 14*11552. 17.77401. 21*22906. 24.58720/ 

DATA <<RJVAL(I. J). J = 1.5). I « 1.9)/ 

1 -0.40276. 0.30012. -0*24970. 0*21836. -0.19647. 

2 0*58187. -0.34613. 0.27330. -0*23330. 0*20701. 

3 0*48650. -0*31353. 0.25474. -0.22088. 0.19794. 

4 0*43439. -0*29116. 0*24074. -0.21097. 0*19042. 

5 0*39965. -0.27438. 0.22959. -0.20276. 0.18403. 

6 0*37409. -0.26109. 0*22039. -0*19580. 0*17849. 

7 0*35414. -0*25017. 0*21261. -0*18978. 0*17363. 

8 0*33793. -0*24096. 0*20588. -0*18449. 0*16929. 

9 0*32438. -0.23303. 0*19998. -0*17979. 0.16539/ 

INPUT PARAMETERS. 

4 READ <5.5000. END ■ 600) (TITLE(I). I « j, 72 > 

READ (5.5001) GAMMA. UE. RLD. ZCOMB. N DROPS. NOZZLE 
IF (GAMMA) 600. 600. 8 
6 READ (5.5004) NJMAX. NONLIN. NEGL. NOUT 
IF (NEGL .EQ. 1) READ (5.5005) SMI. SM2 


L<10). N< 10). NAME(IO). S<10). Sd(lO). TITLE<80). 
RJR00TC10. 5). RJVALOO.5). CK20.20). C(3.20.20). 
0(20.20.20). AMPL(IO). PHASE(IO). AZI (2). 
BESK9.9.9). BES2(9.9.9). BES3<9.9.9). 

V(2). JC( 20). TS( 3. 20). TSQ(20). KMAX(4) 

CRSLT. Cl. ZEJ. ZEPi. ZEF2. CZE. CAZ. CRAD. 

Gl. DCOEF. CGAM. CAX. BOO). BC< 10). YNOZ( 10). 
CNOBMUO). CSSQ(IO). TANINT( 2). RAD1NTO). 

AXINT( 4. 3). CC(4. 10. 10). CD1 ( 10. 10. 10). 

CD2( 10. 10. 10). AX<4). Tl« T2. Dl. D2. D3. D4. 

CD3( 10. 10. 10). CD4<10. 10. 10) 

B /BLK2/ M ( 10). NS(10) 
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c 


c 


IF (NOZZLE •EQ'i 1 ) GO TO 5 
COMPUTE ADMITTANCE FOR QUASI 
Y * (GAMMA - 1.0) * UE/(2.0 
DO 3 J « 1 * NJMAX 


-STEADY NOZZLE. 
* GAMMA) 


AMPL(J) = y 
PHASE(J) = 0.0 
3 CONTINUE 
GO TO 7 


5 DO 6 1 * 1# NJMAX 

READ (5*5003) J* AMFL(J>* PHASE(J) 

6 CONTINUE 

7 DO 10 I * 1* NJMAX 

READ (5*5002) J* L(J), M(J)* N(J)* 
10 CONTINUE 


NS( J)* NAME(J) 


DO 12 J *> 1* NJMAX 

THETA s PHASE(J) * PI/ 180.0 

YR ■ AMPL(J) * COS( THETA) 

YI * AMPL(J) * SIN( THETA) 
YNOZ(J) b CMPLX(YR*YI ) 

12 CONTINUE 
C 


C 

c 

c 

c 


ZE * 2»0 * RLD 
CZE ■ CMFLX(ZE*0.0) 

CGAM » CMPLX( GAMMA* 0.0) 

CAX » CGAM 

IF (NDROPS . EQ. 1 ) CAX « CGAM ♦ 


( 1 * 0 * 0 . 0 ) 


************m**+i'mm m m,n 


ASSIGN ARRAYS FOR ROOTS OF BESSiL 
DO 20 J = 1, NJMAX 

IF ( (M( J) .EO. 0) .AND. (N(J) .EQ. 
MM = M( J) +1 ~ 

NN b N( J) 

S(J) = RJROOT(MM*NN) 

SJ(J) b RJVAL(MM*NN) 


GO TO 25 
15 S(J) b 0.0 
SJ(J) B 1.0 
25 SSQ = S(J) * S( J) 

CSSQ(J) b CMPLX(SSQ*0.0) 
20 CONTINUE 


FLMCTIONS. 
0) ) GO TO 


15 


CALCULATE AXIAL ACOUSTIC EIGENVALUES. 

MAXIMIW VALUES OF L(J), M(J), and N(J) 


c 

c 


c 

c 

c 

c 

c 


LMAX - 0 
MM AX - o 
NMAX > 0 

DO 30 J b i, NJMAX 
IF CL<J) .GT. LMAX) LMAX = L<J) 

IF (M<J) .GT. MMAX) MM AX = M(J) 

IF <N<J) .GT. NMAX) NMAX «= N(J) 

IF CN<J) .NE. N< 1 ) ) KN * 1 

30 CONTINUE 

LMAX «= LMAX + 1 

MMAX ■ MMAX + 1 


COMPUTE EIGENVALUES. 

DO 40 J * L NJMAX 
LL ■ L( J) 

SMN * S<J) 

YAMPL « AMPLCJ) 

YPHASE = PHASE( J) 

B<J) m 1 CRSLT LL ' SMN ' GAMMA,ZE ' YAMPL ' Yp HASE, CRSLT) 

BCCJ) * CON JG ( CRSLT) 

40 CONTINUE 


CALCULATE LINEAR COEFFICIENTS* 


C 

C 


C 

C 


C 

C 


DO 100 NJ * 1* NJMAX 
DO 100 NP ■ J, NJMAX 


105 


ZERO COEFFICIENT ARRAYS. 
DO 105 KC = 1 , 4 
CC(KCsNJ*NP) n C 0*0# 0*0) 
CONTINUE 


112 


ORTHOGONALITY PROPERTY OF TANGENTI A 
IF C NS(NP) .NE. NS(NJ) ) GO TO 10 
IF (M(NP) .NE. M<NJ) ) GO TO 100 
IF (M(NJ) «EQ« 0) GO TO 112 
AZ ■ PI 


GO TO 120 

IF ( NS(NJ) .EQ. |) GO TO 100 
AZ « 2.Q * PI 


EIGENFUNCTIONS. 


120 

125 


ORTHOGONALITY PROPERTY OF 
IF (N<NP) .NE. N(NJ) ) GO 
IF (SCNP)) 125# 122# 125 
SOM • M(NJ) * MCNJ) 

SSO ■ SCNP) * SCNP) 

SJSQ ■ SJ(NJ) * SJ(NJ) 


RADIAL 
TO 100 


EIGENFUNCTIONS. 
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RAD * C SSQ - SQM) * SJSQ/C 2.0 * SSQ) 

GO TO 127 
122 RAD « 0.5 

CALCULATE AXIAL INTEGRALS. 

127 DO 130 NOPT = 1, 4 

CALL AXIAL 1 (NOPT* NP.NJj UE*Z E* Z COMB* CRSLT) 

AXCNOPT) = CRSLT 
130 CONTINUE 

EVALUATE FUNCTIONS AT NOZZLE END. 

ZEJ = CC0SH( CI*BCCNJ)*CZE) 

ZEP1 = CCOSHC CI*BCNP)*CZE) 

ZEP2 » Cl * BCNP) * CSINHC CI*BCNP)*CZE) 

CAZ * CMPLXC AZ*0«0) 

CRAD = CMPLX C RAD* 0*0) 

COEFFICIENT OF THE SECOND DERIVATIVE OF A(P>. 

CCC 1*NJ*NP> = AXC 1 ) * CAZ * CRAD 

COEFFICIENT OF ACP). 

CC(2*NJ*NP> « CCSSQCNP)*AXC 1) - AXC 2) + ZEP2*ZEJ) * CAZ * CRAD 

COEFFICIENT OF THE FIRST DERIVATIVE OF ACP). 

CCC 3*NJ*NP) = C CAX+AXC 3) ♦ C2.0* 0»0)*AXC4) 

1 + CGAM*YN0ZCNP)*ZEP1*ZEJ) * CAZ * CRAD 

COEFFICIENT OF THE RETARDED DERIVATIVE OF ACP). 

CCC 4*NJ*NP> * CGAM * AXC 3) * CAZ * CRAD 

100 CONTINUE 

NORMALIZE LINEAR COEFFICIENTS. 

DO 140 NJ ■ 1* NJMAX 
CNORMCNJ) « CCC 1*NJ*NJ) 

DO 140 NP = 1* NJMAX 
DO 140 KC = 1* 4 

CCCKC*NJ*NP) = CCCKC*NJ*NP) /CNORMCNJ) 

140 CONTINUE 


********************* ******************************************* 

COMPUTE NONLINEAR COEFFICIENTS. 

IF CNONLIN .EQ. 0) GO TO 402 

G1 = CCGAM - Cl. 0*0.0)) * C0»5#0.0) 

COMPUTATIONS OF BESSEL INTEGRALS WHEN ALL SERIES TERMS HAVE THE 
SAME RADIAL MODE NUMBER NCJ). 
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IF (KN .EQ. I) GO TO 170 
DO 150 MP « 1# MM AX 
DO 150 HQ « 1* MM AX 
DO 150 MO * 1. MM AX 
BES l(MP.MQ.MO) > 0*0 
BES2(MP.MQ.M0) * 0*0 
BES3CMP.MQ.M0) « 0*0 
LI « MP - 1 
L2 « MO - 1 
L3 ■ M0 - 1 
LM - LI ♦ L2 
LN - LI ♦ L3 
MN « L2 + L3 

IF ( (L3.EQ.LM) .OR. (L2.EQ.LN) .OR. (Ll.EQ.MN) ) GO TO 160 
GO TO 150 

160 IF (NMAX .EQ. 0) GO TO 165 
A1 » ROROO T(MP. NMAX ) 

A2 « ROROO T(MQ. NMAX) 

A3 ■ ROROO TCMO. NMAX) 

GO TO 167 
165 A1 - 0.0 
A2 - 0.0 
A3 - 0.0 

167 CALL RADIALO.L1.L2.L3.A1.A2. A3. RESULT) 

BESKMP.MQ.MO) - RESULT 

CALL RADI AL( 2.L 1.L2.L3.A1.A2. A3. RESULT) 

BES2(MP.MQ.M0) ■ RESULT 

CALL RADI ALC3.L1.L2.L3.A1.A2. A3. RESULT) 

BES3(MP.MQ«M0) » RESULT 
150 CONTINUE 

170 DO 200 NO - 1. NOMAX 
DO 200 NP * 1. NOMAX 
DO 200 NQ « 1* NOMAX 

CDKNO.NP.NQ) « <0. 0.0*0) 

CD2(N0.NP.NQ) « (0. 0.0*0) 

DO 210 0 - 1. 2 

CALL AZIMTLCO.NP.NQ.NO. RESULT) 

AZI (0) ■ RESULT 

TANINT(O) « CMPLXC RESULT. 0.0) 

210 CONTINUE 

IF (AZI(1>) 220. 225. 220 
225 IF CAZIC2)) 220. 200. 220 

220 IF (KN .EQ. 0) GO TO 222 
LI « M(NP) 

L2 - M (NQ) 
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L3 ■ MCNJ) 

A1 « SCNP) 

A2 » SCNQ) 

A3 • SCNJ) 

GO TO 244 
C 

222 MP ■ M(NP) + 1 
MO • MCNQ) ♦ 1 
MJ ■ MCNJ) + 1 

RADINTC 1) = CMPLXCBES1 CMP,MQ,MJ),0«0) 
RADINTC2) = CMFLXCBES2<MP,MQ,MJ),0.0) 
RADINTC3) = CKPLXCEES3CMP,MQ,MJ),0.0) 

C 

244 DO 240 J * 1, 3 

IF CKN .EQ. 0) GO TO 242 

CALL RADIAL C J,L1,L2,L3,A1, A2, A3, RESULT) 

RADINTC J) = CMPLXC RESULT, 0*0) 

242 DO 240 NC = 1,4 

CALL AXIAL2C J,NC,NP,NQ,NJ,ZE, CRSLT) 
AXINTCNC, J) = CRSLT 
240 CONTINUE 


DO 250 J « 1,4 

T1 ■ G1 * CSSO(NP) * AXINTCJ, 1) 

T2 = G1 * AXINTCJ,3) 

D1 » AXINTCJ, 1) * TANINTC1) * RADINTC3) 

D2 * AXINTCJ, 1 ) * TANINTC 2) * RADINK2) 

D3 * AXINTCJ, 2) * TAN INK 1) * RADINTC 1 ) 

D4 « (T2 - Tl) * TANINTC1) * RADINTC 1 ) 

DCOEF » <0*5,0. 0) * <D1 + D2 + D3 ♦ D4)/CN0RMCNJ) 

IF CJ »EQ. 1) CD1CNJ,NP,NQ) = (1.0>*1>0) * DCOEF 

IF CJ .EQ. 2) CD2CNJ,NP,NQ) a Cl. 0,1.0) * DCOEF 

IF CJ .EQ. 3) CD3CNJ,NP,NQ) = C 1*0, 1*0) * DCOEF 

IF CJ .EQ. 4) CD4CNJ,NP,NQ) a <1.0, -1*0) * DCOEF 

250 CONTINUE 
200 CONTINUE 

**************************************************************** 

CALCULATE COEFFICIENTS FOR EQUIVALENT REAL SYSTEM. 

402 DO 350 NJ ■ 1, NJMAX 
NEWJ = <2 * NJ) - 1 
NEVJ1 a NEVJ ♦ 1 
DO 350 NP ■ 1, NJMAX 
NEWP * <2 * NP) - i 
NEWPl a NEWP ♦ 1 

COEFFICIENTS OF LINEAR TERMS. 
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CCR ■ REAL C CCC 1,NJ,NP) ) 

CCI « AIMAGCCCC 1,NJ,NP) > 
C1<NEWJ,NEWP> - CCR 
C1<NEWJ,NEWP1) ■ -CCI 
Cl (NEWJ1,NEWP> - CCI 
C1<NEWU1,NEWP1 ) ■ CCR 
DO 360 KC = 1,3 
CCR « REAL< CCCKC+ 1,NJ,NP) ) 
CCI » AIMAG ( CCCKC+ 1,N«J,NP) ) 
CCKC,NEWU,NEWP) « CCR 
CCKC,NEWU,NEWP1 ) ■ —CCI 
CCKC,NEWJ1,NEWP) - CCI 
CCKC»NEWU1,NEWP1) ■ CCR 
360 CONTINUE 


C 

C 


COEFFICIENTS OF NONLINEAR TERMS* 
IF CNONLIN *EQ* 0) 60 TO 350 
DO 370 NO ■ 1", NJMAX 
NEWQ » <2 * NO) - 1 
NEW 01 ■ NEWQ ♦ l 


COIR 
CD II 
CD2R 
CDSI 
CD3R 
CD3I 
CD4R 
CD4I 


REAL C CD 1 (NJ,NP,NQ> > 

AIMAGCCD1 <NJ,NP,NQ) ) 

REALCCD2CNJ,NP,NQ) ) 

AIMAGCCD2CNJ,NP,NQ)> 

REAL < CD3C NJ, NP,NQ> ) 

AIMAGCCD3<NJ,NP,NQ>> 

REAL < CD4 C N J, NP, NO)) 

AIM AG ( CD4 ( N J, NP, NO > ) 

DCNEWU,NEWP,NEWQ) » COIR ♦ CD2R ♦ CD3R ♦ CD4R 
D(NEWJ,NEWP,NEWQ1 ) - -CD1I + CD8I - CDsI * tLi 

DCNEWJ,NEWP1,NEVQ> « -coil - C08I ♦ CD 3 I + CDat 

D(NEWJ,NEVP1,NEVQ1 ) = -CD1R ♦ CMR ♦ SsR 
D(NEWJ1,NEWP,NEWQ) * CD1I + CD8I + CD3I ♦ CD4I 
D(NEWJ1,NEWP,NEVQ1) = COIR - CD2R + CD3R - cn/iR 

D<NEWU1,NEWP1,NEWQ) « CD1R + CD2R - CD3R - CD4R 

DCNEW«J1,NEWP1,NEWQ1 ) ■ -CD1I ♦ CD2I ♦ CD3I - rn^r 
370 CONTINUE CD31 CD41 

350 CONTINUE 




COMPUTE COEFFICIENTS FOR THE EQUATION 
IN THE SECOND DERIVATIVES* 


WHICH ARE DECOUPLED 


DO 405 KC « 1 , 4 
KM AX CKC) = o 
405 CONTINUE 


CALCULATE INVERSE OF THE MATRIX 
UMAX “ NoMAX 


C1CI,J>. 
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NJMAX « 2 * NJMAX 


V<1) = l 

CALL GJR( C 1* 20* 20*N JMAX* 0* S 500* JC* V) 

USE INVERSE TO CALCULATE DECOUPLED COEFFICIENTS* 

DO 410 NP = 1 * NJMAX 

LINEAR COEFFICIENTS. 

DO 420 NJ b l* NJMAX 
DO 420 KC c l, 3 
TS(KC*NJ) » 0*0 
DO 420 K = 1 * NJMAX 

TS(KC*NJ) » TS(KC*N J) ♦ CUNJ/K) * C(KC*K*NP) 

420 CONTINUE 

DO 430 NJ = 1 , NJMAX 
DO 430 KC b i, 3 
C(KC*NJ*NP) b TS(KC*NJ) 

ABSVAL b ABS( C(KC*NJ*NP) ) 

IF (ABSVAL *GE» SMI) KMAX(KC) = KMAX(KC) ♦ 1 
430 CONTINUE 

NONLINEAR COEFFICIENTS* 

IF (NONLIN *EQ. 0) GO TO 410 
DO 415 NO b 1 , NJMAX 
DO 440 NJ b NJMAX 
TSQ(NJ) b o*0 
DO 440 K = x, NJMAX 

TSQ(NJ) b TSQ(NJ) + C1(NJ*K) * D(K*NP*NQ) 

440 CONTINUE 

DO 445 NJ b 1 , NJMAX 
D(NJ*NP*NQ) b TSQ(NJ) 

ABSVAL b ABS ( DC N J* NP* NQ ) ) 

IF (ABSVAL *GE* SM2) KMAX(4) - KMAXC 4 ) ♦ 1 
445 CONTINUE 
415 CONTINUE 

410 CONTINUE 


**************************************************************** 

OUTPUT. 

IF (NOUT «GE« 2) GO TO 455 
PRINTED OUTPUT. 

WRITE (6*6001) (TIILE(I)* I » 1 * 72) 

WRITE (6*6002) GAMMA* UE* RLD* ZCOMB 
IF (NDROPS .EQ. 0) WRITE (6*6020) 
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IF (N DROPS .EQ. 1) VRITE (6*6021) 

IF (NOZZLE .EQ. 0) VRITE (6*6012) 

VRITE (6*6004) 

DO 310 0 = 1* OH AX 

VRITE (6*6003) NAME(O)* 0* L(0)* H(0>* N(J)* NS(O)* 
1 S(0>* SJ(O)* B(0)* YNOZ(O) 

310 CONTINUE 

IF (NONLIN .EQ. 0) VRITE (6*6013) 

C 

C OUTPUT OF LINEAR COEFFICIENTS. 

DO 320 KC « 1* 3 

IF (KC .EQ. 1) VRITE (6*6005) 

IF (KC .EQ. 2) VRITE (6*6006) 

IF (KC .EQ. 3) VRITE (6*6007) 

VRITE (6*6008) (0* O > 1* NOMAX) 

VRITE (6*6014) 

DO 320 NO = 1* NOMAX 

VRITE (6*6009) NO* (C(KC*NO*NP>* NP » 1* NOMAX) 

320 CONTINUE 
C 

C OUTPUT OF NONLINEAR COEFFICIENTS. 

IF (NONLIN .EQ. 0) GO TO 452 
DO 400 NO - 1* NOMAX 
VRITE (6*6010) NO 
VRITE (6*6011) (O* 0 « 1* NOMAX) 

VRITE (6*6015) 

DO 400 NP b 1* NOMAX 

VRITE (6*6009) NP* (D(NO*NP*NQ)* NQ « 1* NOMAX) 


400 

CONTINUE 





452 

IF 

(NOUT .EQ. 

0) 

GO 

TO 

4 

455 

IF 

<NOUT . EQ. 

3) 

GO 

TO 

480 


VRITE COEFFICIENTS ON FASTRAND FILE. 

VRITE (9*7001) GAMMA* UE* ZE* ZCOMB* NDROPS* NOMAX 
DO 450 0 « 1* OMAX 

VRITE (9*7002) 0* L(0)* M(0>* N(0>* NS(0>* S(0)* S0(0)* 

1 NAME(O) 

450 CONTINUE 

DO 457 O b i, OMAX 

VRITE (9*7006) O* YNOZ(O)* B(0) 

457 CONTINUE 

DO 460 KC - 1* 3 
VRITE (9*7003) KMAX(KC) 

DO 460 NO ■ 1, NOMAX 
DO 460 NP b i, NOMitX 
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ABSVAL ■ ABS( C(KC*NJ*NP) ) 

460 CONTINUE^* * GE * WWTE <9 ' 7004> NJ. NP, C(KC,NJ,NP) 

WRITE (9*7003) KMAX(4) 

IF (NONLIN .EQ. 0) GO TO 4 
CO 470 NJ = 1* NJMAX 
DO 470 NP = it NJMAX 
DO 470 NQ = 1, NJMAX 
ABSVAL » ABS(D(NJ*NP,NQ)) 

470 CONTINUE^ * 6E * SM2> WR1TE <9*7005) NJ* NP* NQ* D(NJ*NP*NQ) 
GO TO 4 


480 


482 


484 


486 


PUNCHED CARD OUTPUT* 

PUNCH 7001 GAMMA, UE* ZE* ZCOMB* NDROPS* NJMAX 
DO 482 J » 1* JMAX 

^ PUNCH 7002 J^LJJ), MCJ )* N<J), NSCJ), S<J>, SJ(J), 
CONTINUE 

DO 484 J a 1, JMAX 

PUNCH 7006 J* YNOZCJ)* B(J> 

CONTINUE 

DO 486 KC ■ 1* 3 
PUNCH 7003 KMAX(KC) 

DO 486 NJ « 1, NJMAX 
DO 486 NP a j, NJMAX 
ABSVAL a ABS( C(KC*NJ*NP) ) 

IF C ABSVAL »GE» SMI) PUNCH 7004 
CONTINUE 


NJ* NP* C(KC* NJ*NP) 


C 

C 


PUNCH 7003 KMAXC4) 

IF (NONLIN .EQ. 0) GO TO 4 
DO 488 NJ = 1* NJMAX 
DO 488 NP a i, NJMAX 
DO 488 NQ a i, NJMAX 
ABSVAL a ABS(D(NJ*NP*NQ) ) 

468 MNt“uE 4L *“* **’ PU,iC,, 7005 [JJ - »«' DCNJ.NP.NO> 
GO TO 4 

ERROR EXIT 
500 IF ( JC< 1 ) ) 5io* 5lo* 520 
510 JC(1) a ABS(JC(1)) 

WRITE (6*6017) JC(i) 

GO TO 4 
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520 WRITE C 6# 6018 ) JC<1> 

GO TO 4 
600 CONTINUE 

**************************************************************** 


FORMAT 

5000 FORMAT 

5001 FORMAT 

5002 FORMAT 

5003 FORMAT 

5004 FORMAT 

5005 FORMAT 

6001 FORMAT 

6002 FORMAT 
1 

6003 FORMAT 

6004 FORMAT 
1 

6005 FORMAT 

6006 FORMAT 
1 

6007 FORMAT 
1 

6008 FORMAT 

6009 FORMAT 

6010 FORMAT 
1 

6011 FORMAT 

6012 FORMAT 

6013 FORMAT 

6014 FORMAT 

6015 FORMAT 

6017 FORMAT 

6018 FORMAT 

6020 FORMAT 

6021 FORMAT 

7001 FORMAT 

7002 FORMAT 

7003 FORMAT 

7004 FORMAT 

7005 FORMAT 

7006 FORMAT 
END 


SPECIFICATIONS* 

C72A1) 

<4F10* 0,215) 

< 51 5# IX# A4) 

<15# 2F10*0) 

(415) 

(2F10.0) 

< 1H1# IX* 72A 1//) 

<2X*8HGAMMA * # F5* 2# 5X# 5HUE « #F5* 2* 5X* 6HL/D « #F8*5* 

5X# 6KZC0MG * #F5*2/> 

< 2X* A4# 51 5* 6F10* 5/) 

< 2X////2X# 29HNAME J L M N N S# 7X# 3HSMN#,3X* 

7H JM < SMN ) # 7X* 3HEPS# 7X# 3HETA* 8X* 2HYR* 8X* 2HYI //) 

< 1H 1* 45H DECOUPLED COEFFICIENT OF B<P)J CCI#J#P)///> 

< 1H 1# 44H DECOUPLED COEFFICIENT OF THE DERIVATIVE OF* 

6H B( P) * # 5X#8HC<2# J# P)///) 

< 1H 1* 39H DECOUPLED COEFFICIENT OF THE RETARDED* 

20H DERIVATIVE OF B<P) * # 5X#6HC< 3# J* P)///) 

C7X* 1HP* 18*91 12) 

( 2X//2X# I 3* 3X* 10F 12* 6) 

< 1H 1# 42H DECOUPLED COEFFICIENT OF BCP> * DBCQ)/DT* 

19H IN EQUATION FOR BC* I 2* 1H ) ///) 

< 7X* 1HQ* 18*9112) 

<2X* 19HQUASI- STEADY NOZZLE/) 

< 2X//2X* 24HLINEAR COEFFICIENTS ONLY) 

<4X* 1HJ) 

<4X, 1HF > 

<1H1,31H OVERFLOW DETECTED* LAST ROW « #15) 

( 1H 1# 34H SINGULARITY DETECTED# LAST ROW « *15) 

< 2X# * DROPLET MOMENTUM SOURCE NEGLECTED*/) 

< 2X# * DROPLET MOMENTUM SOURCE INCLUDED*/) 

< 4F 10* 5* 21 5) 

< 51 5* 2F 10 • 5# IX* A4) 

< I 5) 

< 21 5#F 1 5*6) 

<3I5*F15*6) 

<15* 4F10.5) 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 


c 

c 


SUBRO UTI ME El G VALCL* SMN# GAMMA* Z E# YAM PL* YPHASE* RESUL T ) 
COMPLEX RESULT 

COMMON /BLK 1/ GSQ* ABSQ* ALBET# SMN SO 


THIS SUBROUTINE COMPUTES THE COMPLEX AXIAL ACOUSTIC EIGENVALUES 
FOR A CYLINDRICAL CHAMBER VI TH A NOZZLE AND STORES THEM IN 
RESULT • 

THE EIGENVALUES ARE COMPUTED BY MEANS OF NEWTONS METHOD. 

THE INPUT PARAMETERS ARE AS FOLLOWS* 

LIS THE AXIAL MODE NUMBER. 

SMN IS THE DIMENSIONLESS ACOUSTIC FREQUENCY. 

GAMMA IS THE SPECIFIC HEAT RATIO. 

ZE IS THE LENGTH- TO- RADI US RATIO. 

YAMPL IS THE NOZZLE AMPLITUDE FACTOR. 

YPHASE IS THE NOZZLE PHASE SHIFT IN DEGREES. 


*****************************11, ***:(, ***,(, 

PI » 3.1415927 
ERR s O.OOOOOOl 

IF C YAMPL) 5# 60# 5 
CALCULATE CONSTANTS. 

5 PHASE * YPHASE * PI/ 180*0 
ALPHA * YAMPL * COSC PHASE) 

BETA b YAMPL * SIN(PHASE) 

GSQ « GAMMA * GAMMA 

ABSQ ■ (ALPHA * ALPHA) - (BETA * BETA) 

ALBET « ALPHA * BETA 
SMNSQ * SMN * SMN 

ASSIGN INITIAL GUESS FOR EIGENVALUE. 

IF (L .EQ. 0) GO TO 45 
RL - L 

PHI e PI/2.0 ♦ PHASE 
XM = RL * PI /ZE 
A b YAMPL/ZE 
XO = XM + A* CO S( PHI ) 

YO = A* SIN (PHI) 

GO TO 47 

45 PHI = PI/4.0 + 0.5* PHASE 
A = YAMPL * 10.0/ZE 
XO = a * COS(PHI) 

YO ■ A * SIN (PHI ) 

ITERATION USING NEWTONS METHOD FOR A SYSTEM OF TWO EQUATIONS 
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IN TWO UNKNOWNS* 

47 LI a 0 
X • XO 

y * yo 

40 CALL FCNS<X#Y#ZE#F#G#FX#FY#GX#GY) 

IF <L1 . EQ« 40) GO TO 50 
RJFG « (FX * GY) - (GX * FY) 

IF (RJFG) 20# 30# 20 
20 DELTAX = <-F * GY ♦ G * FY)/RJFG 
DELTAY = C-G * FX ♦ F * GX)/RJFG 
LI * LI ♦ I 
X » X ♦ DELTAX 

Y ■ Y ♦ DELTAY 
C 

c TEST FOR CONVERGENCE* 

GO T0 B 10 Da ' TAX> ’ GE * *° R * ABS<DELTAY > »GE. ERR) GO TO 

C 

C WARNING MESSAGES 

30 WRITE (6# 6005) 

GO TO 10 

50 WRITE (6# 6006) 

GO TO 10 
C 

C CASE OF HARD WALL (YAMPL » 0)* 

60 RL » L 

X - RL * PI/ZE 

Y ■ 0*0 
C 


C 

C 


10 RESULT «* CMPLX(X#Y) 


6005 

6006 


FORMAT 

FORMAT 

FORMAT 

RETURN 

END 


SPECIFICATIONS. 

C2X//2X#16H JACOBI AN IS ZERO//) 
< 2X//2X# 35HFAILED TO CONVERGE 


IN 40 ITERATIONS//) 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 


SUBROUTINE FCNS<X,Y,ZE,F,G,FX,FY,GX,GY> 


THIS SUBROUTINE COMPUTES THE FUNCTIONS F<X,Y> 
AND THEIR PARTIAL DERIVATIVES WITH RESPECT TO 


AND G<X,Y> 
X AND Y. 


COMMON /BLK 1/ GSQ* ABSQ* ALBET* SMNSQ 

Sd P ™eir H sqIare“°" E ™ IC ruNCTI0NS * ™ E hyperbolic functions 


I = 1 

ARGX » ZE * X 
ARGY « ZE * Y 
10 SX = SIN (ARGX) 

CX * CO S( ARGX) 

SHY = SINH(ARGY) 

CHY « COSH (ARGY) 

IF (I »EQ. 2) GO TO 20 
SXSQ = SX * SX 
CXSQ « CX * CX 
SHYSQ = SHY * SHY 
CHYSQ = CHY * CHY 
ARGX = 2.0 * ARGX 
ARGY « 2.0 * ARGY 
I = 2 
GO TO 10 


COMPUTE TRANSCENDENTAL FUNCTIONS AND THEIR DERIVATIVES 

20 FF = (SXSQ * CHYSQ) - (CXSQ * SHYSQ) 

GG = (CXSQ * CHYSQ) - (SXSQ * SHYSQ) 

HH = 0.25 * SX * SHY 

FFX = ZE * SX * CHY 

GGY * ZE * CX * SHY 

FFY ■ -GGY 

GGX ■ -FFX 

HHX » 0*5 * GGY 

HHY = 0»5 * FFX 

COMPUTE FACTORS 

XYSQ * (X * X) - (Y * Y) 

XY = X * Y 

SMNXY = SMNSQ + XYSQ 

FI = (ABSQ * SMNXY) - (4-0 * ALBET * XY) 

F2 = (ALBET * SMNXY) + (ABSQ * XY) 

G1 • (ABSQ * SMNXY) ♦ (4*0 * ALBET * XY) 

FX1 *» (2.0 * X * ABSQ) - (4*0 * ALBET * Y) 

FX2 - (2.0 * X * ALBET) (ABSQ * Y) 

FY1 = (-2.0 + Y * A3SQ) - (4*0 * ALBET * X) 

FY2 * (-2.0 * Y ♦ ALBET) ♦ (ABSQ * X) 
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rvJ * , (2 1° * X * ABSQ) + <4.0 * ALBET ♦ Y) 
GY1 - <-2. 0 * Y * A BSC) ♦ <4.0 * ALBET ♦ X) 

COMPUTE F<X*Y) AND G(X«Y) 

F = <XYSQ * FF) - <4.0 * XY * HH) 

1 ♦ GSQ * <<F1 * G6) + <4*0 * F2 * HH)) 

G - <XYSQ * HH) + <XY * FF) 

1 ♦ GSQ * <<F2 * GG) - <G1 * HH)) 

COMPUTE THE PARTIAL DERIVATIVES OF F AND G 


1 

2 

3 


FY 


FX « <2.0 * X * FF) ♦ CXYSQ * FFX) 

-4.0 * <<Y * HH) + <XY * HHX) ) 

♦ GSQ + < <FXI * GG) ♦ <F1 * GGX) 

+ <4.0 * FX2 * HH) + <4*0 * F2 * HHX)) 

■ <-2.0 * Y * FF) ♦ <XYSQ * FFY) 

-4.0 * <<X * HH) + <XY * HHY) ) 

♦ GSQ * < <FYI * GG) + <FI * GGY) 

♦ <4.0 * FY2 * HH) + <4.0 
GX » <2.0 * X * HH) ♦ <XYSQ * HHX) 

♦ <Y * FF) + <XY * FFX) 

♦ GSQ * < < FX2 * GG) + <F2 
-<GX1 * HH) - <G 1 * HHX)) 


1 

2 

3 


1 

2 

3 


* F2 * HHY)) 


* GGX) 


GY - <-2.0 * Y * HH) ♦ <XYSQ * HHY) 
l ♦ <X * FF) ♦ <XY * FFY) 

5 + GSQ * <<FY2 * GG> ♦ <F2 * GGY) 

I -<GYI * HH) - <G1 * HHY)) 

RETURN 

END 


¥ 
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SUBRO UTINE AXIAL1(N0PT#NP#NJ#UE#ZE#ZCQMB# RESULT) 


THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 
<0#ZE) OF THE FOLLOWING FUNCTIONS ACCORDING TO THE VALUE 
OF NOPT* 

NOPT « 1 Z (NP) * ZC(NJ) 

NOPT « 2 ZPP(NP) * ZC(NJ) 

NOPT ■ 3 UP * ZCNP) * ZC(NJ) 

NOPT * 4 U * ZP(NP) * ZC(NJ) 

IN THE ABOVE EQUATIONS: 

Z<NP> IS THE AXIAL ACOUSTIC EIGENFUNCTION OF INDEX NP. 

Z<NJ> IS THE AXIAL ACOUSTIC EIGENFUNCTION OF INDEX NJ. 

ZC IS THE COMPLEX CONJUGATE OF THE AXIAL EIGENFUNCTION. 

ZP AND ZPP ARE THE FIRST AND SECOND DERIVATIVES OF THE 
AXIAL EIGENFUNCTIONS RESPECTIVELY. 

U IS THE STEADY STATE VELOCITY DISTRIBUTION AND UP IS ITS 
AXIAL DERIVATIVE. 

THE VELOCITY DISTRIBUTION IS COMPUTED BY THE SUBROUTINE UBAR. 


REAL MAG 

COMPLEX Cl# CZE# BP# BJ# Tl# T2# CH# FI# F2# F3# CZ# ARG# 
1 SI# S2# S3# RESULT# FUNCTC500)# BOO) 

COMMON B 

Cl * C0*0#1.0> 

CZE ■ CMPLX(ZE#0.0) 

BP « BCNP) 

BJ * CONJG(BCNJ) ) 

IF <NOFT .GT. 2) GO TO 50 

CALCULATE INTEGRALS BY MEANS OF ANALYTICAL EXPRESSIONS FOR 
NOPT = l AND NOPT * 2. 

ARG « (BP ♦ BJ) * Cl 
MAG «• CABS < ARG) 

IF <MAG) 20# 25# 20 
20 Tl - CSINH( ARG*CZE)/ARG 
GO TO 30 
25 Tl « CZE 

30 ARG • <BP - BJ) * Cl 
MAG • CABS < ARG) 

IF CMAG) 35# 40# 35 
35 T2 ° CSINHC ARG*CZE) /ARG 
GO TO 45 
40 T2 = CZE 

45 RESULT * <TI ♦ T2) * (0.5# 0*0) 

IF (NOPT .EQ. 2) RESULT » -B(NP) * BCNP) * RESULT 
GO TO 100 
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C NUMERICAL EVALUATION OF INTEGRALS FOR NOPT = 3 AND NOPT - 4. 
c 

C COMPUTE STEP SIZE FOR SIMPSON INTEGRATION* 

50 N = 50 
RN * N 

RESULT = <0. 0/0*0) 

IC * ZCOMB 
IC « 2 - IC 
C 

DO 90 J » l, IC 

IF CJ *EQ* 1) H - ZCOMB * ZE/RN 

IF CJ *EQ* 2) H ■ <1.0 - ZCOMB) * ZE/RN 

IF (J *EQ* 1) 20 « 0*0 

IF CJ «EQ* 2) ZO » ZCOMB * ZE 

NPI ■ N + 1 

CH « CMPLXCH/O.O) 

C 

C COMPUTE INTEGRANDS. 

DO 60 I * l, NPI 
STEP =1-1 
Z ■ C STEP * H) + ZO 

IF ((I.EQ.l) .AND* CJ.EQ.2)) Z * Z ♦ H/ I 00*0 
IF CNOPT *EQ. 3) CALL UBARC2/UE/ZE/ZC0MB/Z/F) 

IF CNOPT .EG. 4) CALL UBARC 1/ UE/ZE/ZCOMB/Z/ F) 

FI * CMPLXCF/O.O) 

CZ ■ CMPLXCZ.O.O) 

ARG • Cl * BP 

IF CNOPT .EG. 3) F2 « CCOSHCARG*CZ) 

IF CNOPT .EQ. 4) F2 « ARG * CSINHCARG*CZ) 

ARG . Cl * BJ 
F3 « CCOSHC ARG* CZ > 

FUNCTCI) ■ FI * F2 * F3 
60 CONTINUE 
C 


PERFORM SIMPSON INTEGRATION. 

NM1 = N - 1 

51 « FUNCTCI) ♦ FUNCTCNP1 ) 

52 ■ C0*0/0*0) 

53 ■ CO* 0/0.0) 

DO 70 I » 2/ N/ 2 

52 ■ S2 ♦ FUNCTCI) 

70 CONTINUE 

DO 80 I « 3/ NM1/ 2 

53 - S3 ♦ FUNCTCI) 

80 CONTINUE 

RESULT = RESULT ♦ 

I CH * CS1 ♦ C4«0/0«0)*S2 ♦ 

90 CONTINUE 

100 CONTINUE 
RETURN 
END 
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C2.0/0.0)*S3)/C3.0/0.0) 


ooooooooooooo 


SUBROUTINE UBARCNOPT* UE*ZE*ZCOMB*Z*RESULT) 

THIS SUBROUTINE CALCULATES THE STEADY STATE VELOCITY 
DISTRIBUTION FOR UNIFORMLY DISTRIBUTED COMBUSTION COMPLETED AT 
Z ■ ZCOMB * ZE WHERE: 

UE IS THE EXIT MACH NUMBER. 

ZE IS THE DIMENSIONLESS LENGTH. 

Z IS THE AXIAL COORDINATE. 

IF NOPT ■ 1 THE DISTRIBUTION IS CALCULATED. 

IF NOPT » 2 THE DERIVATIVE IS CALCULATED. 

IF NOPT = 3 THE SECOND DERIVATIVE IS CALCULATED. 


ECZ = ZCOMB * ZE 



GO 

TO 

( 10*20*30) 

* NOPT 



10 

IF 

CZ 

• LE. 

ECZ ) 

RESULT 

= 

UE * Z/ECZ 


IF 

CZ 

• GT. 

ECZ) 

RESULT 

ss 

UE 


GO 

TO 

AO 





0 

CM 

IF 

CZ 

•LE. 

ECZ) 

RESULT 

at 

UE/ECZ 


IF 

CZ 

• GT. 

ECZ) 

RESULT 

e 

0.0 


GO TO 40 
30 RESULT = 0.0 
40 CONTINUE 
RETURN 
END 
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SUBROUTINE AZIMTLCNOPT,NP,NQ,N«J, RESULT) 

DIMENSION NFCNC 3 ), SGC2) 

COMMON /BLK2/ MCIO), NSC 10 ) 



THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 

0F ™ E FOLL OWING FUNCTIONS ACCORDING TO THE VALUE 


NOPT * 1 THCNP) * THCNQ) * THCNJ) 

NOPT » 2 THPCNP) * THPCNQ) * THCNJ) 


IN THE ABOVE EQUATIONS* 

THCNQ)* AND THCNJ) ARE THE TANGENTIAL EIGENFUNCTIONS 
AND NP, NQ, AND NJ ARE THEIR INDICES* 

THP IS THE DERIVATIVE OF THE TANGENTIAL EIGENFUNCTIONS* 


IF NS « 1 TH *= SINCM+THETA) 
IF NS - 2 TH « CO SCM* THETA) 



RESULT >0*0 
FACTOR = 1.0 
PI ■ 3*1415927 


C 


C 


DISTINGUISH BETWEEN 
DO 10 K1 - 1* 3 
NFCNCK1) « l 
10 CONTINUE 

IF CNSCNJ) • EQ* 2) 

IF CNOPT *EQ. 2 ) 

IF CNSCNP) *EQ*2) 

IF CNSCNQ) • EQ* 2) 

GO TO 30 

20 IF CNSCNP). EQ.l) 

IF CNSCNQ). EQ.l) 

DO 40 HI » 1,2 
SGCK1) =1*0 
IF CNFCNCK1) 

40 CONTINUE 

FACTOR * SGC 1 ) * 


SINES AND COSINES. 


NFCNC3) * 
GO TO 20 
NFCNC 1 ) = 
NFCNC2) = 

NFCNC 1 ) «= 
NFCNC 2 ) a 


2 

2 

2 

2 


EQ. 1 ) SGCK 1 ) « -1.0 


SGC2) * MCNP) * MCNQ) 


30 N SUM « 0 

DO 50 K 1 = 1 , 3 
NSUM «■ NSUM + NFCNCK 1 ) 
50 CONTINUE 
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IF C (NSIM »EQ. 3) •OR* CNSUM *EQ* 5)) GO TO 60 
IF (NSUM .EQ. 4) GO TO 70 
IF CN SUM .EQ. 6) GO TO 80 

70 KOPT = 2 

IF CNFCNC1) .EQ. 2) GO TO 72 
GO TO 74 
72 LL = MCNP) 

MM = MCNQ) 

NN = MCNJ) 

GO TO 90 

74 IF CNFCNC2) .EQ. 2) GO TO 76 
GO TO 78 
76 LL * M(NQ) 

MM a MCNP) 

NN a MCNJ) 

GO TO 90 
78 LL = MCNJ) 

MM a MCNP) 

NN a MCNQ) 

GO TO 90 

80 KOPT a i 
LL a MCNP) 

MM a MCNQ) 

NN a MCNJ) 

COMPUTE VALUES OF THE INTEGRALS. 

90 IF CCLL.NE.O) .AND. CMM.NE.O) -AND. CNN.NE.O)) GO TO 101 
GO TO 103 

101 LM « LL + MM 

LN a LL + NN 

MN a MM + NN 

IF CCNN.EQ.LM) .OR. CMM.EQ.LN)> RESULT a PI/2.0 
IF CLL .EQ. MN) GO TO 102 
GO TO 104 

102 IF CKOPT .EQ. 1) RESULT a PI/2.0 
IF CKOPT .EQ. 2) RESULT a -PI/2.0 
GO TO 104 

103 IF CCLL.EQ.O) .AND. CMM.EQ.O) .AND. CNN.EQ.O)) GO TO 105 

IF C CKOPT. EQ. 1 ) .AND. CNN.EQ.O) .AND. CLL.EQ.MM>) RESULT ■ PI 
IF CCKOPT.EQ.l) .AND* CMM.EQ.O) .AND. CLL.EQ.NN)) RESULT * PI 
IF C CLL .EQ. 0) .AND. CMM .EQ. NN)) RESULT a pi 
GO TO 104 

105 IF CKOPT .EQ. 1) RESULT a 2.0 * PI 

104 CONTINUE 

RESULT a FACTOR * RESULT 
60 CONTINUE 
RETURN 
END 
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SUBROUTINE RADI AL(NOPT#L#M#N# A# B# C# RESULT) 

THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 
<0#1) OF THE FOLLOWING PRODUCTS OF THREE BESSEL FUNCTIONS* 

NOPT = I JL( A*R> * JM(B*R) * JN(C*R) * R 

NOPT » 2 JLCA*R) * JM(B*R) * JN(C*R)/R 

NOPT ■ 3 JPL(A*R) * JPM(B*R> * JN(C*R) * R 

JL IS THE BESSEL FUNCTION OF FIRST KIND OF ORDER L 

JPL IS THE DERIVATIVE OF JL WITH RESPECT TO R 
L# M# N ARE NON-NEGATIVE INTEGERS 
A# B# C ARE REAL NUMBERS 

DIMENSION FUNCT( 200) 

DOUELE PRECISION DN# DH# DSTEP# DR# ARGI# ARG2# ARG3# 

1 BESI# BES2# BES3# BESH# BESL# PROD# 

2 FUNCT# BESLIM# SI# S2# S3 

NN « 100 
DN * NN 
DH ■ 1 .0/DN 
NPl ■ NN 1 

DO 10 I - I# NPl 
DSTEP = I - 1 
DR ■ DH * DSTEP 
ARGI a A * DR 
ARG2 » B * DR 
ARG3 a c * DR 

CALL <JBES(N# ARG 3# BES3# S500) 

IF (NOPT .EQ. 3) GO TO 101 
CALL «JBES(L# ARG 1# BESI# S500) 

CALL JBES(M#ABG2#BES2#S500) 

GO TO 102 

101 IF <L .EQ. 0) GO TO 103 

CALL JBES CL+ 1 # ARGI# BESH# S 500) 

CALL JBES(L-1# ARGI# BESL# S500) 

BESI « A * (BESL - BESH)/2.0 
GO TO 104 

103 CALL JBES( 1# ARG 1#BES1# $500) 

BESI a -BESI * A 

104 IF (M .EQ. 0) GO TO 105 
CALL JBES(M+ 1# ARG2# BESH# S500) 

CALL JBESCM- 1, ARG2# j<ESL# $500) 

BES2 « B ♦ (BESL - BESH)/2.0 
GO TO 102 


105 CALL UBES( 1# ARG2#BES2# S500) 

BES2 = -BES2 * B 
102 PROD = BES1 * BES2 * BES3 

IF (NOPT .EQ. 2) GO TO 110 
FUNCT(I) = PROD * DR 
GO TO 10 

110 IF <1 .EQ. 1) GO TO 111 
FUNCTtl) = PROD/DR 

GO TO 10 

111 BESLIM = 0.0 

IF ((L.EQ.l) .AND. (M.EQ.O) .AND* (N.EQ.O)) BESLIM ■ A/2.0 

IF ((L.EQ.O) .AND. (M.EQ.l) .AND. (N.EQ.O)) BESLIM = B/2.0 

IF ((L.EQ.O) .AND. (M.EQ.O) .AND. (N.EQ* 1)) BESLIM = C/2.0 

FUNCT(I) « BESLIM 
10 CONTINUE 

NM1 a NN - 1 

51 - FUNCT(l) ♦ FUNCT(NPl) 

52 » 0*0 

53 » 0*0 

DO 20 I « 2# NN# 2 

52 ■ S2 + FUNCT(I) 

20 CONTINUE 

DO 30 I = 3# NM1# 2 

53 » S3 ♦ FUNCT(I) 

30 CONTINUE 

RESULT 8 DH * (SI ♦ 4«0*S2 ♦ 2.0*S3)/3»0 
GO TO 501 

500 WRITE (6# 6000) 

6000 FORMAT ( 1H1# 10HERROR JBES) 

501 CONTINUE 
RETURN 
END 


93 


n o 


i 


C 

c 

c 

c 

c 

c 

c 

c 

c 

e 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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SUBROUTINE AXI AL2CN0PT*NC0N J*NP* NQ*NJ*ZE* RESULT) 


THIS SUBROUTINE CALCULATES THE INTEGRAL OVER THE INTERVAL 
<0*ZE) OF THE FOLLOWING FUNCTIONS ACCORDING TO THE VALUES 
OF NOPT AND NCONJ: 

FOR NCONJ * 1 AND: 

NOPT « 1 ZCNP) ♦ Z (NQ) * ZCCNJ) 

NOPT ■ 2 ZPCNP) * ZPCNQ) * ZCCNJ) 

NOPT « 3 ZPPCNP) ♦ ZCNQ) * ZCCNJ) 

FOR NCONJ «= £ AND: 

NOPT * 1 ZCNP) * ZCCNQ) * ZCCNJ) 

NOPT « 2 ZPCNP) * ZPCCNO) * ZCCNJ) 

NOPT * 3 ZPPCNP) * ZCCNQ) * ZCCNJ) 

FOR NCONJ « 3 AND: 

NOPT e 1 ZCCNP) * ZCNQ) * ZCCNJ) 

NOPT * 2 ZPCCNP) * ZPCNQ) * ZCCNJ) 

NOPT « 3 ZPPCCNP) * ZCNQ) * ZCCNJ) 

FOR NCONJ *» 4 AND: 

NOPT *= 1 ZCCNP) * ZCCNQ) * ZCCNJ) 

NOPT e 2 ZPCCNP) * ZPCCNQ) * ZCCNJ) 

NOPT = 3 ZPPCCNP) * ZCCNQ) * ZCCNJ) 


IN THE ABOVE EQUATIONS: 

ZCNP)* ZCNQ)* AND ZCNJ) ARE THE AXIAL ACOUSTIC 
AND NP* NQ* AND NJ ARE THEIR INDICES. 


EIGENFUNCTIONS 


ZP IS THE FIRST DERIVATIVE OF THE AXIAL EIGENFUNCTIONS* 

ZPP IS THE SECOND DERIVATIVE OF THE AXIAL EIGENFUNCTIONS. 

ZC AND ZPC ARE COMPLEX CONJUGATES OF Z AND ZP RESPECTIVELY. 


REAL MAG 

COMPLEX Cl* CF* CZE* BP* BQ* BJ* SIM* RESULT* 
1 ARG C A ) * FUNCTC 4)* BCIO) 

COMMON B 


CALCULATE INTEGRALS BY MEANS OF ANALYTICAL EXPRESSIONS. 
Cl ■ CO. 0*1.0) 


CF = <0.25*0.0) 

CZE « CM PLX C Z E* 0 . 0 ) 

BP * BCNP) 

BQ = BCNQ) 

BJ ® CONJGCBCN J) ) 

IF CCNCONJ .EQ. 2) .OR. CNCONJ .EQ. 4)) 
IF CNCONJ *GT. 2) BP « CONJGCBP) 

ARG Cl) « CBP ♦ BQ + BJ) * Cl 


BQ » CONJGCBQ) 
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12 

15 

10 


30 

50 


ARG ( 2 > = (BP + BQ - BJ) * Cl 
ARGO) * (BP - BQ BJ) * Cl 
ARG ( A ) • CBP - BQ - BJ) * Cl 
DO 10 J = 1,4 


MAG * CABSCARG(J)) 

IF (MAG) 12, 15j 12 

FUNCT(J) = CSINH(ARG(J)*CZE)/ARG(J) 

n n th t a 


GO TO 10 
FUNCT(J) 
CONTINUE 
IF (NOPT 


CZE 


•EQ. 2) 


GO TO 30 

SUM = FUNCT(l) + FUNCT( 2) + FUNCTO) 
RESULT ® CF * SUM 

IF (NOPT »EQ« 3) RESULT = -BP * BP * 
GO TO 50 


SUM = FUNCT(l) + FUNCT( 2) - 

RESULT e -cf * BP * BQ * SUM 

CONTINUE 

RETURN 

END 


FUNCTO) 


+ FUNCT( 4) 
RESULT 
- FUNCTO) 


APPENDIX D 


PROGRAM LCYC3D: A USER'S MANUAL 


General Description 

Using the three-dimensional second-order theory described in this report 
Program LCYC3D calculates the nonlinear stability characteristics of a cylin- 
drical combustion chamber with distributed combustion and a conventional noz- 
zle. The response of the burning rate to pressure oscillations is described 
by Crocco's time-lag model. For given values of the operating parameters 
(i.e.> n 5 T, y, u^, and L/d), a given series expansion, and a given initial 
disturbance Program LCYC3D integrates Eqs. (C-38) to obtain the time behavi- 
or of the unknown mode -amplitude functions (i.e., B.(t)). From this infor- 
mation a time history of the pressure oscillation is determined. The program 
determines the final amplitude of the pressure oscillation attained in a 
linearly unstable engine (i.e., limit-cycle anplitude) . Since the second- 
order analysis does not predict "triggering", however, the threshold ampli- 
tude above which a finite amplitude disturbance can trigger instability in a 
linearly stable engine (i.e., triggering limit) is not calculated by Program 
LCYC3D. For either transient or limit-cycle conditions, the program prints 
out time histories of both pressure and axial velocity perturbations from 
which the amplitude, frequency, and wave shapes can be deter mi ned The op- 
-i° n to produce plotted output using a CALCOMP plotter is also provided. 

Program Structure 

A flow chart for Program LCYC3D is given in Fig. (D-l) . This program 
performs the following operations: (l) reads the input data, (2) calculates 
the initial conditions, ( 3 ) numerically integrates the differential equations, 

(4) tests for limit cycles (optional), and ( 5 ) prints and plots the resulting 
solutions. 

The inputs to the program include the data generated by Program C0EFFS3D, 
the combustion parameters n and t, various control numbers, and a description 
of the initial disturbance. The data from C0EFFS3D is read first and then 
printed out. Next the space dependent coefficients appearing in the series 
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Figure D-l. Flow Chart for Program LCYC3D 












expansions for § Q , and are computed and printed out. These coeffi- 
cients are calculated by Subroutine PHICFS for use in the computation of the 
pressure and axial velocity perturbations. The remaining input data is then 
read, and following program execution, control is returned to this point (see 
Fig. D-l) so that several cases (i.e., different values of n and t) may be 
run for a given set of coefficients generated by C0EFFS3D. 

After input of the initial amplitudes of the real parts (i.e., Bg . -^(t)) 
of the complex amplitude functions, the initial amplitudes of the imaginary 
parts (i.e., are calculated such that the nozzle admittance condition 

is satisfied for -T s t s; 0. These amplitudes are then printed out. Next 
the integration step-size, At, is calculated such that the interval 
-T s t s 0 is divided into NDIV equal increments. Assuming a sinusoidal 
initial disturbance, the initial amplitudes of B 2J _ 1 (t) and B (t) are used 
to calculate these functions and their derivatives at each of the NDIV + 1 
discrete points in -T ^ t ^ 0. These values are needed in order to start the 
numerical solution of the differential equations (i.e., Eqs. (C-38)). The 
initial values of the amplitude functions are stored in the array U(l,j) 
where the index I varies from 1 (t = -T) to NDIV + 1 (t = 0) and the index 
J identifies the function. The corresponding initial values of the pressure 
and velocity perturbations are then printed out. This section also calcu- 
lates the coefficients C 2 (j,p) - nC^(j,p) and nC^(j,p) which are the coeffi- 
cients of dB^/dt and d^B p (t - r)J/dt in Eqs. (C-38). 

After the starting values are calculated, Eqs. (C-38) are solved using 
a modified form of the fourth order Runge Kutta method. Starting at t = 0 
(I = NlJlV+l) , the amplitude functions at t + At are calculated, using the 
Subroutine RHS to evaluate the functions . . .B^) on the right hand 

sides of Eqs. (C-38). The amplitude functions and the coefficients from 
PHICFS are then used to compute the pressure and axial velocity perturbations 
by Subroutine PRSVEL. The values of the amplitude functions at t + At are 
stored in U(l + 1,J), while the pressure and axial velocity perturbations are 
stored in the arrays PRESS(NPRES) and AXVEL(NPRES) where NPRES specifies the 
locations in the chamber where the data is calculated. Pressure data at one 
location (specified by NLOC) is also stored in the array PRS(l + l) . After 
checking for maximum and minimum values of U(l,j) and PRS(l), the data may 
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be printed out (if NTEST = 0 and TSTART <: t <; TQUIT) or stored in plot arrays 
as desired. The time is then increased by At (i.e., I is increased by l) and 
the calculations are repeated. This process continues until 250 integration 
steps have been computed (t = 250 At) , after which transfer is made to the 
limit-cycle section. 

In the limit-cycle section a test for a limit-cycle is made if NTEST = 1. 
If the test is satisfied, NTEST is set to zero so that no further tests will 
be made and the results can be printed or plotted. In either case the final 
values (for 250-NDIV £ I ^ 250) replace the initial values (for 1 ^ I ^ 
NDIV+l) in the arrays U(l,j) and PRS(l), I is again assigned the value NDIV+1, 
and another 250 integration steps are calculated. This process continues un- 
til one of the following conditions is satisfied: (l) NTEST = 0 and t > TQUIT, 

(2) a limit-cycle is reached and t > TQUIT, and (3) more than 250 cycles of 
the pressure oscillation have been computed (MAXNO > 500) . At this point the 
numerical calculations are terminated and the time history of the pressure 
amplitude (maxima and minima) are printed out and/or plotted as desired. 

As can be seen from Fig. D-l the output is not confined to a single 

section of the program but is produced in several different sections. Thus 

data is printed out or plotted shortly after it is calculated, which greatly 

reduces the amount of core storage required. All plots are generated by Sub- 

20 

routine GRAPHS which uses standard Univac 1108 plot routines. 

FORTRAN listings of Program LCYC3D and Subroutines PHICFS, PRSVEL, RHS, 
and GRAPHS are provided at the end of this appendix. 

Input Data 

A precise definition of the input data required to run the computer 
program is given below. This input data consists of three parts: (l) the 

control number N0UTCF, (2) the parameters and coefficients generated by Pro- 
gram C0EFFS3D and (3) the data describing the cases to be run (see Fig. D-l). 
For each input case the following information must be provided: (l) the 

combustion parameters n and t; (2) a series of control numbers; and (3) infor- 
mation describing the initial disturbance. 

The control number N0UTCF determines whether the coefficients from 
C0EFFS3D will be printed, and it appears on the first card of input. This 
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card is followed by the coefficient deck generated by C0EFFS3D and the data 
describing the cases to be run. Since the coefficient data has already been 
described in Appendix C, it will be omitted from the following detailed des- 
cription of the input. As in Appendix C the location number refers to the 
columns of the card. Again three formats are used for input: "A" indicates 

alphanumeric characters, "i" indicates integers, and "F" indicates real num- 
bers with a decimal point. For the "I" formats the values are placed in 
fields of five locations, while a field of ten locations is used with the 
"F" formats. In either case the numbers must be placed in the rightmost 
locations of the allocated field. 

No. of 


Cards 

Location 

Type 

Input Item 

Comments 

i 

1-5 

I 

NOUTCF 

If 0: coefficients are not 


printed out. 

If 1: linear coefficients 
only are printed out . 

If 2: all coefficients are 
printed out . 


1-72 

A 

TITLE 

Title used to label plots. 

1-10 

F 

EN 

Interaction index, n. 

11-20 

F 

TAU 

Time-lag, t. 

21-30 

F 

H 

Time-increment for numerical 
integration, ^t. * 

31-40 

F 

T START 

Time at which output of solu- 
tions begins. 

41-50 

F 

TQUIT 

Time at which output of solu- 
tions ends. 

1-5 

I 

NTEST 

If 0: compute transient beha- 
vior . 




If 1: compute limit-cycle be- 
havior. 

6-10 

I 

JM0DE 

Identifies the amplitude func 


tion used to test for limit- 
cycles . 


* This value is adjusted slightly by the program to divide the interval 
~T s t s 0 into NDIV equal parts. 
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No. of 



Cards 

Location 

Type 

Input Item 

Comments 



11-15 

I 

NLOC 

Determines location for wall 
pressure maxima and minima. 

If 1: z = 0, e = 0° 

If 2: z = 0, e = 45° 






If 3: z = 0, e = 90° 



16-20 

I 

NTERMS 

Number of amplitude functions 
given initial values. 



21-25 

I 

NPZ 

Determines how secondary 
instability zones are handled. 






If 0: all instability zones 

retained . 






If 1: secondary zones elim- 

inated . 



26-30 

I 

NOUT 

Determines output. 






If 0: printed output only. 






If 1 s NOUT ^ 6: both print- 

ed and plotted output, NOUT 
gives number of last plot 
produced . 


If Is 

NOUT <; 6 the 

following two cards are 

read: 


1 

1-10 

F 

YHl(l) 

Maximum ordinate for pressure 
plots . 



11-20 

F 

YHI(5) 

Maximum ordinate for velocity 
plots . 



21-30 

F 

YLAB(l) 

Interval for ordinate labeling 
of pressure plots. 

* 


31-40 

F 

YLAB(5) 

Interval for ordinate labeling 
of velocity plots. 

* 

1 

1-5 

I 

ITICY(l) 

Number of ordinate tic marks 
for pressure plots. 



6-10 

I 

ITICY(5) 

Number of ordinate tic marks 
for velocity plots. 



11-15 

I 

NFIRST 

Gives the number of the first 
plot produced. 



16-20 

I 

NOMIT 

If 0: amplitude plot produced 

If 1: amplitude plot omitted. 
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No. of 
Cards 

Location Type 

Input Item 

Comments 


End of 

input for 1 £ NOUT £ 6. 




NTERMS 

1-5 I 

j 

Identifies complex amplitude 
function. 



6-15 I 

AST 

Amplitude of sin(uut) term in 
initial conditions. 

*■ 


16-25 I 

ACT 

Amplitude of cos (cut) term in 
initial conditions. 



The input data describing the cases to be run is given on a series of 
three or more cards. These cards are preceded by a title card which gives 
a title (TITLE) to be used to identify any plots produced by the run. This 
title appears before the first plot generated and does not appear on the 
printed output. The title card is included only for the first case of the 
run; on all subsequent cases it is omitted. 

The first card of the series gives the interaction index, n, and the 

time-lag, T, for the motor under consideration (EN and TAU); the time- 

increment, At, used in the numerical integrations (H); and the times (TSTART 

and TQUIT) at which output begins and ends. For all cases considered in this 

report a time -increment (dimensionless) of H = 0.050 was used, which gives 

about 70 steps per cycle for the IT mode. For T = 1.7 this input value was 

adjusted by the program to obtain H = 0.04857 which divides -T £ t s 0 into 

35 equal parts. For transient cases (NTEST = 0) printed output is given for 

TSTART s t s TQUIT. When the limit-cycle behavior is calculated (NTEST = l), 

TSTART and TQUIT are measured from the time at which the limit-cycle is 

reached, t T _. Thus the limit-cycle solutions are printed out for 

(t Tri + TSTART) £ t £ (t + TQUIT). Two or three cycles of limit-cycle data 
LC LC 

for the IT mode are obtained with TSTART = 0 and TQUIT = 10. For plotted 
output, the time axis is always 10 units long, therefore (TQUIT - TSTART) > 10 
to obtain plots. 

The second card of the series gives the control numbers, NTEST, JMODE, 
NLOC, NTERMS, NPZ, and NOUT. The task to be performed by Program LCYC3D is 
specified by NTEST. If NTEST = 0 the transient behavior (growth or decay) of 
the pressure oscillation is determined, while for NTEST = 1 the program 
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searches for a limit-cycle amplitude. JMODE identifies the "principal" 
series term, the amplitude function used in the limit-cycle test. This is 
usually the lowest frequency mode (i.e., IT or 1L) in the approximating 
series expansion. NLOC gives the location at which the amplitude- time his- 
tory (maxima and minima) of the wall pressure perturbation is calculated. 

The number of complex series terms A (t) receiving initial values is speci- 
fied by NTEEMS, while all other series terms are initially zero. The para- 
meter NPZ determines how the secondary instability zones (phantom zones) are 
handled by Program LCYC3D. For NPZ = 1 the phantom zones are eliminated by 

dropping the combustion terms for a given mode when t > T where 1 

cut 


T = = 2TT 

cut W 


[ 


? ,2 2 

S 2 

mn 2 


] 


1 

2 


(D-l) 


A similar procedure was used in the axial instability studies by Lores and 
Zinn. The transverse instability data presented herein was obtained with 
NPZ = 0, while NPZ = 1 was used in the axial instability studies to facili- 
tate comparison with the results of Ref. (3). The last control marker NOUT 
determines which plots, if any, are produced. For NOUT = 0 no plots are 
produced. For 1 £ NOUT ^ 6, NOUT gives the number of the last plot produced, 
where the plots are numbered as given in Table D-l below: 


Table D-l. Numbering of Plots. 


No. of Plot 
( NPLOT ) 

Quantity 

Plotted 

Axial 

Location 

Azimuthal 

Coordinate 

1 

Pressure 

Injector 

0° 

2 

IT 

IT 

45° 

3 

11 

II 

VO 

o 

o 

4 

II 

Nozzle 

0° 

5 

Axial Velocity 

II 

0° 

6 

Nozzle Boundary 
Term 

It 

0° 
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The nozzle boundary term given on the last plot is discussed- later in this 
appendix. 

If plots are produced, two additional cards are needed to give the max- 
imum and minimum values of the variables to be plotted, YHl(NPLOT) and 
YLO(NPLOT); the intervals for ordinate labeling ( YLAB ( NPLOT ) ) ; and the num- 
ber of ordinate tic marks, ITICY(NPLOT) . All of the plots are symmetric >» 

about the time-axis so that YLO(NPLOT) = -YHl(NPLOT), and ITICY( NPLOT) must 
be negative to obtain the centerline. Since the ordinate scales and labeling * 
are the same for all pressure plots (NPLOT = 1,2, 3, 4) this data is read for 
NPLOT = 1 only; likewise the data for the last two plots is read for NPLOT = 

5 only. In addition NFIRST gives the number of the first plot produced, 
giving additional control over the number of plots produced. NOMIT deter- 
mines whether a plot of pressure amplitude versus time (location specified 
by NLOC) is produced. 

The remaining cards give the initial amplitudes of the complex series 
terms, A..(t), needed to start the numerical integration. Only the anplitudes 
of the real parts, B^^-^t), are S-'- ven on these cards, while the amplitudes 
of the imaginary parts, B^(t), are determined from the nozzle admittance 
condition. For each value of J the anplitudes AST and ACT are assigned to 
the arrays AS(NP) and AC(NP) where NP = 2J - 1. The computation of the 
amplitudes of the imaginary parts, AS(NP + l) and AC(NP + l) , is discussed 
later. The initial values of the series terms are then calculated from the 
formula: 

B (t) = AS(NP)sin(u) t) + AC(NP)cos(oj t) (-T £ t £ 0) (D-2) 

where Ul p is the acoustic frequency. The derivatives, dB /dt, are also re- 
quired for starting the numerical integration; they are obtained simply by • 
differentiating Eq. (D-2) . 

The proper input for pure standing and pure spinning single-mode initial 
disturbances is given as follows. For a standing mode, only the cos(me) terms 
are retained in the series and NTERMS = 1. A single card is read giving the 
amplitude of the initial disturbance. For a spinning mode, both sin(m0) and 
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cos(mQ) terms are included in the series expansion. It is convenient to pair 
these terms such that the index J corresponds to a sin(me) term and J + 1 
corresponds to a cos(m0) term. For an initial disturbance of amplitude A 
spinning in the counterclockwise direction (9 increasing), NTERMS = 2 and 
two cards are read giving the following data: 


J : AST = A and ACT = 0 

(D-3) 

J + 1 : AST = 0 and ACT = A 


In both cases above initial amplitudes are required only for the mode initi- 
ally present, and the initial amplitudes of all other modes included in the 
series expansion are zero. 

The proper input for Program LCYC3D will be illustrated with the follow- 
ing example. Assuming that the velocity potential § is expressed in terms of 

•X* 

the 1R, IT, and 2T modes , it is desired to determine the limit-cycle beha- 
vior of a linearly unstable engine (n = 0.57486, T = 1.7? u g = 0.2, L/D = 0 . 5 ) 
with a nozzle admittance of A = 0.02 and cp = 4-5°. Sample input is given for 
the case of a spinning IT mode disturbance of amplitude 0.3. The principal 
series term is the cos(m0) term for the IT mode (i.e., BQ^(t)), thus JM0DE 
= 2. Plots are desired for the pressure, axial velocity, and nozzle boundary 
condition at the nozzle entrance, thus N0UT = 6 and NFIRST = 4. 

To run the case described above the data deck must be assembled as 
follows. The card specifying N0UTCF is followed by the coefficient deck 
produced by Program C0EFFS3D; in this example it contains the information 
given in the sample output for C0EFFS3D shown in Appendix C. The coefficient 
deck is followed by the data for the case to be run as shown in the sample 
input below: 


This is the same case used to illustrate Program C0EFFS3P. 
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Table E>-2. Sample Input. 


lii 

II3II 

1 

U 

1 

MU 

II 

HI 

1 

■XI 

1 

mu 

1 

MU 

1 

IU 

1 

KJ 

1 

■a 

I 

EJ 

1 

■a 

1 

EJ 

I 

El 

1 

EH 

1 

0 

1 

0 

1 

□ 

1 

□ 

1 

0 

1 

□ 

1 

ED 

1 

0 

1 

0 

1 

□ 

1 

□ 

1 

a 

1 

□ 

1 

a 

1 

B 

■ 

n 

■ 

0 

fl 

B 

| 

n 

| 

B 

fl 

B 

| 

B 

fl 

b 

| 

B 

1 

Kt^FI 

|| 

B 

| 

H 

n 

mm 

Ml 

B 

B 

B 



n 

ns 

B 

B 

B 

B 

n 

B 


B 


E3 

Sf 

10 


i 

i 

IQ 

1 

i 

OS 

D 

a 

S 

Q 

G 

fl 

1 

0 

B 

B 

■ 

B 

3 

1 

D 

OB 

B 

B 

a 

1 

a 

a 

a 

0 

a 

a 

5 

fl 

fl 

fl 

fl 

| 

I 


1 


rr 

n 

n 

n 

B 

B 

B 

B 

B 



B 

B 

K3 

IS 

n 

19 

B 

B 

E 

B 

B 

B 

E 

m 

1 

■1 


1 

IS 

IB 

1C 

n 

la 

II 

II 

II 

I 

II 

■ 

■ 

IQ 

II 

B 

II 

II 

II 

II 

II 

II 

IG 

II 

IG 

a 

II 

■ 

■ 

■ 

■ 

■ 

II 

IQ 

II 

IE 

|| 

|| 

| 

|| 

|| 

■ 

IQ 

0 

II 

B 

a 

n 

wm 

SI 

B 

B 


ES 

B 

B 

n 

B 

B 

B 

B 

B 

B 

B 

E 


sS 

1 

11 


0 

1 

■ 

■ 

■ 

a 

■ 

■ 

1 

1 

0 

1 

1 

1 

1 

B 

1 

1 

1 

1 

G 

1 

■ 

1 

■ 

0 

1 

1 

1 

1 

1 

1 

1 

■ 

■ 

1 

■ 

| 

| 

fl 

| 

1 

I 

I 

fl 

H 

n 

n 


B 


B 

B 

B 

B 

B 

B 

B 

E 



1 


1 

IIElflEjfll 

■ 

II 

1 

G 

1 

G 

U 

1 

1 

1 

1 

1 

1 

1 


a 

1 

1 

1 

1 

1 

1 

B 

fl 

B 

a 

I 

1 

1 

I 

| 

1 

1 

1 

| 

1 

c 


Q 

B 

B 

B 

B 

Fl 

M 

HIE 

\mw 

_ 

a 

a 

0 

■ 

■ 

■ 

1 

u 

1 

1 

1 

1 

E 

1 

1 

1 

■ 

1 

1 

1 

■ 

1 

1 

1 

1 

1 

1 

1 

1 

■ 

■ 

■ 

1 

■ 

1 

fl 

fl 

| 

fl 

fl 

fl 

fl 

fl 

a 

KM 

SI 

n 

n 

D 

B 


■Ml 

m 

m 

m 

Fl 

Pi 

PI 

m 

Fl 

— 

II 


ll 

s 





■ 

■ 

1 

0 

1 

a 

1 

■ 

1 

1 

1 

1 

1 

G 

1 

G 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

■ 

1 

1 

1 

■ 


n 

Oi 


11 

a 

a 

fl 

n 


Fl 

rnmmm 

11 


III 



■ 






El 

II 

01 


■ 

■ 

■ 

1 

1 


El 


SI 

■ 

■ 


3 

j 

n 

j 

3 

3 

3 

: 

3 

3 

3 

3 

3 

3 

3 

j 

: 

: 

n 

: 

3 



Coefficients in Series for 


*t> V ^ V 


As seen from Eq_. (13) the real parts of the time and space derivatives 
of the velocity potential (i.e., 9 , , 9 r ? 9 5 9 ) are needed in order to com- 

U X 0 2 

pute the pressure perturbation. Differentiating the complex series expan- 
sion given by Eq. ( 9 ) and evaluating at the chamber wall (r = l) gives the 
following expansions: 


N jn 

r-> dA 

It * 

P=1 


z p (z)o p Ce)R p (i) 


N 

£ c t (v,z,Q) 
P=1 


dA 

dt 


$ 


e 


N 


I y^y^p^Hpd) 

p=i 


N 

£ c 0 (p,z,e)A p (t) 
p=i 


(d-4) 


(D-5) 


(D-6) 


N N 

$ = y A (t)z'(z)® (e)R (1) = y C (p,z,e)A (t) 
zZjPPPP p 

P=1 P=1 


where the complex coefficients C, , C , and C are functions of z and e. The 

"C 0 z 

quantity, $ r , is not needed since = 0 at the chamber wall. The complex 

coefficients Cl , C , and C are calculated by Subroutine PHICFS and are 
t 0 z 

assigned to the variables, Cl, C2, and C3 respectively. The coefficients in 
the series expansions for the corresponding real parts (i.e.,9^, cp Q , 9 z ) are 
related to the complex coefficients by: 


C t (2p-l, z, e) = Re^C t (p, z, e)J 

C t (2p, z, 9) = -Im[^C t (p, z, e)J 


(D-7) 


/ / 

where similar relations hold for C Q and C z - The real coefficients are stored 
in the arrays CFT(NPRES, NP) , CFTH(NPRES, NP) , and CFZ(NPRES, NP) where NPRES 
determines the location in the chamber as given in Table D-3 below: 


Table D-3. Chamber Locations for Pressure Calculations. 


Axial 

NPKES Location (z) 


Azimuthal 
Location (q) 


10 0 ° 

2 0 45° 


3 0 90° 

4 z e 0° 

5 z e 45° 

6 z 90° 

e 


Initial Amplitudes 

The initial amplitudes of the real parts of the complex series terms 
(i.e., B (t)) are specified in the input to the program. The initial 
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amplitudes of the imaginary parts (i.e., B 2J (t)), however, are calculated 
such that the nozzle admittance condition is satisfied for -T ^ t sO. This 
is done by introducing the linear expressions for u' and p' into the nozzle 
admittance relation and assuming periodic solutions. This yields a set of 
linear algebraic equations relating the amplitudes of the real and imaginary 
parts of the complex series terms. For given values of the amplitudes of the 
real parts, AS(NP) and AC(NP), these equations are solved to obtain the anpli- 
tudes of the imaginary parts, AS(NP + l) and AC(NP + l) . The following for- 
mulas are used in this calculation. 


where 


AS(NJ + 1) = -(r 2&1 - r x a 2 ) / (a^ + a|) 
AC(NJ + 1) = (r i&1 + r 2 a 2 ) / (a^ + a|) 


r 1 = a 2 [ac(HJ)] - a 4 [as(NJ)] 
r 2 = " a 4 [ AC ( NJ )] - a 3 [as(NJ)] 


(D-8) 


(D-9) 


and 


a l - (l + yY r U e )CFZ(NPRES, NJ+l) - yY . tu . CFT ( NPRES , NJ+l) 

1 J 

a 2 = ^^(NPHSS, NJ+l) + yY u CFZ( NPRES, NJ+l) 

(D-10) 

a 3 - "(! + y Y r u e )CFZ( NPRES, Nj) + yY.OU .CFT (NPRES, NJ) 

X J 

a 4 = yY r «M. CFT (NPRES, Nj) + yYuCFZ( NPRES, Nj) 

In Eqs. (D-8) through (D-10) is the acoustic frequency and CFT and CFZ are 
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the coefficients in the series for cp and cp computed previously. The above 
conditions are applied at a pressure anti-node for each series term, there- 
fore NPRES = 4 (z = z g , 0 = 0°) for a cos(me) term and NPRES = 6 (z = z , 0 = 90°) 
for a sin(m0) term. 

For nozzles with phase shifts of cp = 90° and cp = 270° the quantity 
2 2 

a l + a 2 vanishes and Eqs . (D-8) become indeterminate. In these cases the amp- 
litudes of the imaginary parts are given by: 


AS(NJ + l) = AC(NJ) 
AC(NJ + 1) = AS(NJ) 


(D-ll) 


which provides a good approximation to the nozzle admittance condition. 

Integration of the Differential Equations 

For purposes of numerical integration Eqs. (C-38) are written as an 
equivalent system of first order differential equations as follows: 



dfh 



/ 



(D-12) 

(D-13) 


where the dependent variables are now B. and B.. These equations are solved 

0 0 

numerically using the fourth order Runge-Kutta method. Due to the presence 
of retarded variables in Eqs. (D-12) and (D-13) the formulas (see Ref. 21 ) 
used in the Runge-Kutta method must be slightly modified. 

The appropriate formulas for applying the Runge-Kutta method to problems 
involving a time-delay are readily obtained by considering a single equation 
of the following form: 


— = f (x,t) + g[x( t - T)] 


(D-14) 
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Noting that at any step of the integration the value of x(t - t) has already 
been determined from previous steps, the function g can be considered to be a 
known function of time g(t). 

Since x(t) is computed only at discrete points x n (t ) it is desired that 
the retarded variable x(t n - t) will coincide with such previously computed 
points. This can be accomplished by choosing the step-size At such that it 
divides the time-lag T into k equal increments. Thus j = k^t and the Runge- 
Kutta formulas which apply to Eq. (D-l4) can now be written as: 

Vl = x n + Kb + 2k 2 + ^3 + b) 
b = { f(x n>V + g(x n-k ) } 4t 

k 2 = { f(x n * k l/ 2 ’ \ + At/ 2 ) + 8 < x „. k+ j)}i t (D-15) 

k 3 = {b x n + V -J - ' t -„ + A'-/ 2 ) + e(x n _ k+ i)}it 

k 4 = { f(x n + k 3 - b + At) + S^n-k+l^A* 

Equations (D-15) are readily extended to handle the system of equations given 
by Eqs. (D-12) and (D-13). It is seen from Eqs. (D- 15 ) that k values of the 
dependent variables prior to the initial values are needed to start the inte- 
gration. 

Although the initial wave shape can be an arbitrary function of time, it 
is assumed that initially the mode -amplitudes are sinusoidal functions of time 
oscillating with the natural frequency u>.. Thus each mode- amplitude function 
is expressed in the following form: 

B . ( t) = AS( J)sin(d) .t) + AC(j)cos(i».t) 

<-) J J 

(D-16) 
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B. 

J 


(t) = w^AS( j)cos(uAt) - AC(j)sin(uut) J 


where -f s t £ 0. 

/■ 

In Program LCYC3D both the functions B.(t) and the derivatives B.(t) are 

stored in the same array U(l,j). The B.(t) (N functions) are stored in the 

J 

first half of the array (l s J sN), while the remaining space (N + 1 ^ J ^ 2N) 

is used to store the values of B^'(t) . Thus for a given value of j (l s j s N), 

B.(t) is stored in U(l,j) and B((t) is stored in U(l,J + N) . In addition the 
0 _ J 

retarded variables B'(t - t) are stored in the array RV(j,K) as follows: 

0 

rv(j,i) = s'(t - ;) 

J 


RV(J,2) = RV(J,3) = B'(t - t + At/2) (D-17) 

J 

rv(j, 4) = B'(t - ; + At) 

o 


' r _ f 

The values of B.(t - T + ^t/2) are computed from B.(t - x )> B.(t - x + At) , 
f _ 0 J J 

and B . ( t - T + 2 At) using a three-point interpolation. 

I 

Pressure and Axial Velocity Perturbations 

From the calculated time dependence of the series terms Program LCYC3D 
computes the dimensionless pressure perturbation, p', with the aid of Eqs. 
(D-4) through (D-6) and either Eq. (13) for NDROPS = 0 or Eq. (A- 6 ) for 
NDROPS = 1. The pressure is calculated at the injector face (z = 0) and the 
nozzle entrance plane (z = z ) for three angular positions along the peri- 
phery of the chamber (i.e., r = 1; 0 = 0 °, 45°, 90°). The results are stored 
in the array PRESS(NPRES) where NPRES gives the location according to Table 
D-3 . The axial velocity perturbation at the nozzle entrance, u', is calcula- 
ted for 0 = 0°, 45°, 90° using the relation u' = cp and Eq. (D-6), and the 
results are stored in AXVEL(K), where K = NPRES-3. In addition the quantity, 
Re£-yY§^J, is calculated at the nozzle entrance for 0=0° and assigned to 
the variable YPHI. From Eq. (2) it is seen that YPHI is the axial velocity 
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at the nozzle entrance (i.e., u') if the nozzle admittance condition is exactly 
satisfied. Since the solutions generated by Program LCYC3D are approximate, 
the difference between u^ and YPHI is a measure of the accuracy of this approx- 
imation at the nozzle boundary. 

Maximum and Minimum Values 

In order to determine the transient behavior and limit-cycle amplitudes 
it is necessary to follow the growth or decay of the amplitudes of the series 
terms and the pressure perturbation. The maxima and minima of the principal 
series term (specified by JMODE) are assigned to the array UMAX(MAXNO) where 
MAXNO is a counter variable. For the pressure perturbation, m a x i mum and mini- 
mum values at the location specified by NLOC are stored in PMAX(MAXP), and the 
corresponding times of maximum and minimum are stored in TIMAX(MAXP). Since 
the solutions are calculated only at discrete points, the maximum and minimum 
values are computed using a three-point interpolation scheme. 

Calculation of Limit-Cycle Amplitude 

A limit-cycle amplitude is calculated by specifying an initial disturbance 
and continuing the step-by-step integration of Eqs. (D-12) and (D-l;j) until a 
periodic solution is obtained; that is, the amplitude of the oscillation 
remains essentially constant. The test for convergence to a limit cycle is 
performed upon a single series term, usually the most important term in the 
series, in the following manner. After the first 500 integration steps, 
usually about 10 cycles for the IT mode, the amplitude of the principal series 
term A.^ is compared with its amplitude after 250 integration steps Aq. If 
the change in amplitude | A^ - Aq | is greater than the maximum permissible 
change e , the calculations are continued and the change in amplitude during 
the next 250 integration steps is calculated. The process is repeated until ' 

I ^k " A k-1 I < e which point the computation is terminated. The amplitudes 
used in the above calculations are determined by averaging the absolute values 
of UMAX(MAXNO) over the last two complete cycles for each 250 integration 

steps. A value of c = 0.001 is used in Program LCYC3D which gives sufficient 
accuracy for most cases. 
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Output 


Printed Output . The printed output produced by Program LCYC3D consists 
of the five sections discussed below. 

Section 1 is a restatement of the input from Program C0EFFS3D. It in- 
cludes the following information: (a) the ratio of specific heats (GAMMA), 

the mean flow Mach number at the nozzle entrance (UE), the dimensionless 
chamber length (ZE), the length of the combustion zone as a fraction of the 
chamber length (ZCOMB), and the number of series terms (real) NJMAX; (b) a 
statement regarding the presence or absence of the droplet momentum source; 

(c) the parameters which describe and identify each term in the series expan- 
sion; (d) the nozzle admittance (YR and Yl) and the axial acoustic eigenvalue 
(EPS and ETA) for each series term; (e) the nonzero linear coefficients, 

C(KC, NJ, NP) ; and (f) the nonzero nonlinear coefficients, D(NJ, NP, NQ) . 

The nonlinear coefficients are omitted from the output for NOUTCF = 1, and 
no coefficients are printed out for NOUTCF = 0. 

Section 2 gives the coefficients needed for computation of the wall 
pressure waveforms; that is, the coefficients in the series for cp^, cp^, and 
cp z - These are given for each of the NJMAX series terms at each of the six 
locations specified by NPKES (see Table D-3). 

Section 3 gives the initial amplitudes (AS(j) and AC(j)) of all series 
terms included in the assumed initial disturbance. This section also states 
whether the limit-cycle behavior is calculated and whether plots are produced. 

Section 4 gives the time-dependent solutions for the following quantities: 
(a) the injector pressure perturbation at 0=0°, 45°, 90 °; (b) the nozzle 
pressure perturbation at 0 = 0°, 45°, 9°°J (c) the nozzle axial velocity 
perturbation at 0 = 0 , 45 , 90 ° and (d) the nozzle boundary term, Re£— yY$^J, 
at 0 = 0°. This output is given in two parts: (l) the initial values for 
-T ^ t £ 0 and (2) the solutions for £ t £ t f , where t^ and t f are deter- 

mined by the input parameters TSTART and TQUIT (see discussion on Input). On 
the first page of each part a heading gives the interaction index, n, and the 
time-lag, and the chamber parameters, y, u g , and L/D. 

Section 5 gives the time history of the pressure amplitude (maximum and 
minimum values) for the chamber location specified by NLOC. This information 
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is printed as an array of number pairs giving the value of the pressure maxi- 
mum or minimum (upper number) and the corresponding time of maximum or minimum 
(lower number). This information is useful in determining the growth (or decay) 
rate of the transient solutions, and it provides a check on the convergence of 
the solution to a limit-cycle. 

Plotted Output . According to the values of NOUT and NFIRST the pressure 
and axial velocity waveforms given in Section 4 of the printed output may be 
plotted using a Calcorap plotter. The data over the dimensionless time interval 
for printed output, tj s t s t f , is plotted in sections of 10 units in length 
beginning at t = t^. Thus for each quantity plotted, N plots are produced 
where N is the largest multiple of 10 contained in the interval t ± £ t s: t . 

The data left over (i.e., for t^ + ION ^ t ^ t^) is not plotted. All quantities 
to be plotted for a given time interval are plotted before proceeding to the 
next time interval. 

The data given in Section 5 of the printed output (pressure maxima only) 
is also plotted if NOUT > 0 and NOMT = 0. The abscissa and ordinate ranges 
for this plot are not specified in the input, but are calculated such that 
all of the data falls within these ranges. This plot is always the last plot 
produced . 

All of the above plots are scaled to fit on standard 8|" x 11" paper 
and scissor-lines are plotted for trimming plots to this size. The data is 
plotted as individual points using a small circle symbol, and all of the values 
computed during the given time interval are plotted. Before the first plot is 
produced the identifying title (see Input) is printed. 

Sample Output . The following sample output illustrates the printed and 
plotted output produced by Program LCYC3D for the sample input given in Table 

D-2. 
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OAMiiA = 1.200 


Table D-4. 


UE = .200 


Sample Output, Section 1. 


ZE = 1.00000 ZCOMB = 1.00 


droplet momentum source is neglected 


NAME 


J L M N NS 


SMN JM (SMN) 


AOll 

eon 

A021 

B021 

B001 


1 

0 

1 

1 

1 

1.84118 

.58187 

z 

0 

1 

1 

2 

1.84118 

.58187 

3 

0 

2 

1 

1 

3.05424 

.48650 

4 

0 

2 

1 

2 

3.05424 

.48650 

5 

0 

0 

1 

2 

3.83171 

-.40276 


J 

YR 

TI EPS 

ETA 

1 

.01414 

.01414 .08122 

.19451 

2 

.01414 

.01414 .08122 

.19451 

3 

.01414 

.01414 .10617 

.25115 

4 

.01414 

.01414 .10617 

.25115 

5 

.01414 

.01414 . H993 

.28170 

NUMBER 

OF COEFFICIENTS Ctl.NJ.NP) IS 

10 

C(l» 1. 

1) = 

3.39060 


Ctl. Zt 

2) = 

3.39060 


Ctl. 3. 

3) = 

3.39060 


Ctl. 4 » 

4) = 

3.39060 


Ctl. 5. 

51 = 

9.33021 


Ctl. 6. 

b) = 

9.33021 


Ctl. 7. 

7) = 

9.33021 


Ctl. 8. 

8) = 

9.33021 


Ctl. 9. 

9) = 14.68491 


Ctl. 10.10) = 14.68491 


number 

OF COEFFICIENTS Ct2.NJ»NPj IS 

10 

C 1 2. 1. 

1) = 

.26153 


C t 2 . 2. 

2) = 

.26153 


Ct2. 3. 

3) = 

.26153 


C 12. 4. 

4) = 

.26153 


C t 2» 5. 

5) = 

.26457 


C ( 2 » 6 . 

o ) - 

•2b457 


C 1 2. 7. 

7) = 

.26457 


Ct2. 8. 

8) = 

.26457 


C 1 2. 9. 

9) = 

.26654 


C 1 2 . 1 0 » 

10) = 

.26654 


number 

OF COEFFICIENTS C 1 3 » N J » NP ) IS 

10 

Cl3. 1, 

1) = 

.24000 


C t 3 . Zt 

2) = 

.24000 


C l 3. 3. 

3) = 

.24000 


C t 3. 4. 

4) = 

.24000 


C t 3. 5. 

5) = 

.24000 


C l 3 . 6. 

b) = 

.24000 


C t 3. 7. 

7) = 

.24000 



NJMaX 


Table D-4. (Continued) 


C(3* 8* 8) = .24000 
C (3* 9* 9> = .24000 
C(3*l0*l0> = .24000 

HUMBER OF COEFFICIENTS D(NJ*NP.NO) IS 


D( 

1' 

if 

7) 

z 

-1.73504 

Dl 

if 

It 

9) 

z 

-2.33866 

D( 

if 

if 

b) 

z 

1.73504 

0( 

if 

5f 

3) 

z 

1.49783 

Dl 

if 

7 r 

i) 

z 

-1.49783 

Dl 

If 

9f 

1) 

z 

-1.96281 

Oi 

2f 

2f 

6) 

z 

-1.73505 

0 i 

2f 

2f 

LO) 

z 

-2.33867 

0( 

2f 

4f 

6) 

z 

1.73505 

Dl 

2f 

6f 

4) 

z 

1.49784 

Oi 

2f 

8f 

2) 

z 

-1.49784 

Oi 

2f 

iOf 

2) 

z 

-1.96282 

Oi 

3f 

lr 

5) 

z 

1,73504 

Oi 

3f 

3f 

7) 

z 

1.73504 

Oi 

3f 

3f 

9) 

z 

-2.33866 

Oi 

if 

5f 

1) 

z 

1.49783 

Oi 

if 

7 f 

3) 

z 

1.49783 

Oi 

3f 

9r 

3) 

z 

-1.9b281 

Oi 

**f 

2f 

6) 

z 

1.73505 

Oi 

Hf 

4 f 

6) 

z 

1.73505 

Oi 

Hf 

4 f 

10) 

z 

-2.33867 

Oi 

4f 

6f 

2) 

z 

1.49784 

Oi 

4 f 

8 f 

4) 

z 

1.49784 

Oi 

4f 

lOr 

4) 

z 

-1.96282 

Oi 

5f 

If 

3) 

z 

-1.13133 

D< 

5 f 

3f 

1) 

z 

-1.13133 

Oi 

5f 

b f 

9) 

z 

-3.07465 

Oi 

5f 

9f 

5) 

z 

-2.81865 

D( 

6 1 

2f 

4) 

z 

-1.13132 

Oi 

6f 

4 f 

2) 

z 

-1.13132 

Oi 

6f 

6f 

iO) 

z 

-3.07469 

Oi 

6f 

iOf 

6) 

z 

-2.81868 

D( 

7f 

if 

1) 

z 

1.13133 

Oi 

7 f 

3f 

3) 

z 

-1.13133 

Oi 

7 f 

7 f 

9) 

z 

-3.07465 

Oi 

7 f 

9 f 

7) 

z 

-2.81865 

Oi 

8f 

2f 

2) 

z 

1.13132 

Oi 

8f 

4 f 

4) 

z 

-1.13132 

Oi 

8 9 

6f 

10) 

z 

-3.07469 

Oi 

6f 

iOf 

8) 

z 

-2.81868 

Oi 

9f 

if 

1) 

z 

1.040B7 

Oi 

9f 

3f 

3) 

z 

1.04087 

Oi 

9f 

5f 

5) 

z 

-.21090 

Oi 

9f 

7f 

7) 

z 

-.21090 

Oi 

9f 

9f 

9) 

z 

4.18784 

Oi 

XOf 

2 f 

2) 

z 

1.04087 

0(10* 

4 f 

4) 

z 

1.04007 

0(10* 

6f 

to) 

z 

-.21091 

0(10* 

8f 

tt) 

z 

-.21091 

0(10* 

IOf 

1U) 

= 

4.18793 


Table D-5. Sample Output, Section 2. 


COEFFICIENTS FOK COMPUTATION 0^ WALL PRESSURE WAVEFORMS 


1 

2 

3 

4 

5 

6 
7 
a 

9 

10 

1 

2 

3 

4 

5 

6 
7 

a 

9 

10 

1 

2 

3 

4 

5 

6 
7 
6 
9 

10 

1 

2 

3 

4 
b 
6 
7 
6 
9 

10 

1 

2 

3 


COEFFICIENTS IN SERIES FOR: 


z 

THETA 
t DEGREES) 

TIME 

DERIVATIVE 

theta 

DERIVATIVE 

AXIAL 

DERIVATIVE 


.000 

.0 

.0000000 

.5818700 

.OOOOOOO 

.000 

• 0 

.0000000 

.0000000 

•OOOOOOO 

.000 

.0 

.5818700 

• ooooooo 

.OOOOOOO 

.000 

• 0 

.0000000 

.0000000 

.OOOOOOO 

.000 

• 0 

.0000000 

.9730900 

•OOOOOOO 

.000 

• 0 

.0000000 

.0000000 

.OOOOOOO 

• 000 

.0 

.4865000 

.0000000 

.OOOOOOO 

.000 

• 0 

.0000000 

.0000000 

•ooooooo 

.000 

• 0 

-.4027600 

.0000000 

.ooooooo 

.000 

• 0 

.0000000 

.0000000 

.0090000 

.000 

45.0 

.4114442 

.4114442 

.ooooooo 

.000 

45.0 

.0000000 

.0000000 

.ooooooo 

.000 

45.0 

.4114442 

-.4114442 

.ooooooo 

.000 

45.0 

.0000000 

.0000000 

.ooooooo 

.000 

45.0 

.4865000 

-.0000000 

.ooooooo 

.000 

45.0 

.0000000 

.0000000 

.ooooooo 

.000 

*5.0 

-.0000000 

-.9730000 

.ooooooo 

.000 

45.0 

.0000000 

.0000000 

.ooooooo 

.000 

45.0 

-.4027600 

.0000000 

.ooooooo 

.000 

45.0 

.0000000 

.000000O 

.ooooooo 

.000 

90.0 

.5818700 

-.0000000 

.ooooooo 

.000 

90.0 

.0000000 

.0090000 

.ooooooo 

.000 

90.0 

-.0000000 

-.5818700 

.ooooooo 

.000 

90.0 

.0000000 

.0000000 

.ooooooo 

.000 

90.0 

-.0000000 

-.9730000 

.ooooooo 

.000 

90.0 

.0000000 

.0000000 

.ooooooo 

.000 

90.0 

-.4865000 

.0000001 

.ooooooo 

.000 

90.0 

• 0000000 

.0000000 

.ooooooo 

.000 

90.0 

-.4027600 

.OOOOOOn 

.ooooooo 

.000 

90.0 

.0000000 

.0000000 

.ooooooo 

1.000 

.0 

.0000000 

.5909578 

.ooooooo 

1.000 

.0 

.0000000 

.0092403 

. ooooooo 

1.000 

.0 

.5909575 

.0000000 

.0181736 

1.000 

• 0 

.0092403 

.OOOOOOO 

.0185766 

1.000 

• 0 

.0000000 

.9981959 

.ooooooo 

1.000 

.0 

.0000000 

.0261690 

.ooooooo 

1.000 

• 0 

.4990979 

.OOOOOOn 

.0251885 

1.000 

.0 

.0130845 

.OOOOOOO 

.0263938 

1-000 

• 0 

-.4158379 

.oooooon 

-^0261428 

1.000 

• 0 

-.0137546 

.ooooooo 

-.0278051 

1.000 

45.0 

.4178700 

.4178700 

.0128507 

1.000 

45.0 

.0065339 

.0065339 

.0131356 

1.000 

45.0 

.4178700 

-.4178700 

.0128507 


vDCfc^jacn-Pt^tv)*— c cd ^ (rn ^ 


Table D-5 . (Continued) 


1.000 

45.0 

.0065339 

-.0065339 

.0131356 

1.000 

45.0 

.4990979 

-.0000000 

.0251885 

1.000 

45.0 

.0130845 

-.nooooon 

.0263938 

1.000 

45.0 

-.0000000 

-.9981959 

-.0000000 

1.000 

45.0 

-.0000000 

-.0261690 

-.0000000 

1.000 

45.0 

-.4158379 

.0000000 

-.0261428 

1.000 

45.0 

-.0137546 

.0000000 

-.0278051 

1.000 

90.0 

.5909575 

-.0000000 

.0181736 

1.000 

90.0 

.0092403 

-.0000000 

.0185766 

1.000 

90.0 

-.OOOOUUO 

-.5909575 

-.0000000 

1.000 

90.0 

-.0000000 

-.0092903 

-.0000000 

1.000 

90. 0 

-.0000000 

-.9981959 

-.0000000 

1.000 

90. 0 

-.0000000 

-.0261690 

-.oooooou 

1.000 

90.0 

-.4990979 

.0000001 

-.0251885 

1.000 

90.0 

-.0130845 

.0000000 

-.0263938 

1.000 

90.0 

-.4158379 

.0000000 

-.0261428 

1.000 

90.0 

-.0137546 

.0000000 

-.0278051 


Table D-6. Sample Output, Section 3- 


INITIAL conditions are of the porm: 

utl,0) = AC(J)*COS«FREQ*T) + AS(J)*SIN(FREO*T))» * EXP(DAMP*T) 


J 

DAMPING 

FREQUENCY 

ACIJ) 

AS(J> 

1 

.00000000 

1.84118000 

.ooonoooo 

.30000000 

2 

.00000000 

1.84118000 

-.30278619 

-.00209447 

3 

.00000000 

1.84118000 

.30000000 

.00000000 

4 

.00000000 

1.84118000 

-.00209447 

.30278619 


THE LIMIT-CYCLE BEHAVIOR IS CALCULATED. 


this run produces plotted output 
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Figure D-3. Sample Amplitude Plot 
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FORTRAN Listing 


*************** PROGRAM LCYC3D ******************************* 

THIS PROGRAM CALCULATES THE NONLINEAR BEHAVIOR OF 
TRANSVERSE* AXIAL* OR COMBINED LONGITUDINAL-TRANSVERSE 
INSTABILITIES IN A CYLINDRICAL COMBUSTION CHAMBER WITH 
UNIFORM PROPELLANT INJECTION* DISTRIBUTED COMBUSTION 
PROCESS* AND A CONVENTIONAL NOZZLE. THE COMBUSTION PROCESS 
IS DESCRIBED BY CROCCO’S TIME-LAG MODEL. BOTH TRANSIENT 
AND LIMIT-CYCLE SOLUTIONS ARE CALCULATED. 

THE FOLLOWING INPUTS ARE REQUIRED: 

(1) THE CONTROL NUMBER* NOUTCF. 

(2) THE COEFFICIENTS FROM PROGRAM COEFFS3D. 

( 3) THE DATA DECK. 

NOUTCF DETERMINES PRINTOUT OF COEFFICIENTS. 

IF NOUTCF = 0 COEFFICIENTS ARE NOT PRINTED OUT* 

IF NOUTCF = 1 LINEAR COEFFICIENTS ONLY ARE PRINTED OUT. 

IF NOUTCF = 2 ALL COEFFICIENTS ARE PRINTED OUT. 

THE DATA DECK CONSISTS OF THE FOLLOWING CARDS: 

FIRST CARD: 

EN IS THE INTERACTION INDEX. 

TAU IS THE TIME LAG. 

HIS THE INTEGRATION STEP SIZE. 

TSTART IS THE TIME AT WHICH OUTPUT STARTS. 

TQUIT IS THE TIME AT WHICH COMPUTATIONS ARE TERMINATED. 

SECOND CARD: 

NTEST IS TASK CONTROL NUMBER: 

IF NTEST = 0 COMPUTE TRANSIENT BEHAVIOR. 

IF NTEST = 1 COMPUTE THE LIMIT-CYCLE BEHAVIOR. 

JMODE IS THE MODE-AMPLITUDE USED TO TEST FOR LIMIT-CYCLES. 

NLOC DETERMINES THE LOCATION OF THE WALL PRESSURE MAXIMA 

AND minima: 

IF NLOC = 1 LOCATION IS Z c 0* THETA = 0 DEGREES. 

IF NLOC = 2 LOCATION IS Z = 0* THETA = 45 DEGREES. 

IF NLOC = 3 LOCATION IS Z = 0* THETA «* 90 DEGREES. * 

N TERMS IS THE NUMBER OF TERMS GIVEN INITIAL VALUES. 

NPZ DETERMINES HOW SECONDARY STABILITY ZONES (PHANTOM 

ZONES) ARE HANDLED. 4 

IF NPZ = 0 PHANTOM ZONES ARE RETAINED. 

IF NPZ = 1 PHANTOM ZONES ARE ELIMINATED. 

NOUT IS THE OUTPUT CONTROL NUMBER. 

IF NOUT = 0 PRINTED OUTPUT ONLY. 

IF NOUT > 0 BOTH PRINTED AND FLOTTED OUTPUT* NOUT 
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DETERMINES THE NUMBER OF THE LAST PLOT 
PRODUCED. 

DATA FOR SETTING UP PLOTS (THIRD AND FOURTH CARDS): 

YHI ( 1 ) IS THE MAXIMUM ORDINATE FOR PRESSURE PLOTS. 

YHI ( 5) IS THE MAXIMUM ORDINATE FOR VELOCITY PLOTS. 

NOTE: THE ORDINATE SCALES FOR PRESSURE AND VELOCITY PLOTS 
ARE SYMMETRIC ABOUT ZERO. 

YLAB IS THE INTERVAL FOR ORDINATE LABELING FOR ABOVE PLOTS. 
ITICY IS THE NUMBER OF ORDINATE TIC MARKS FOR ABOVE PLOTS* 
NOTE: ITICY SHOULD BE NEGATIVE FOR PRESSURE AND VELOCITY FLOTS 
TO OBTAIN CENTERLINE. 

NFIRST IS THE NUMBER OF THE FIRST PLOT PRODUCED. 

NOMIT DETERMINES WHETHER AMPLITUDE FLOT IS PRODUCED: 

IF NOMIT = 0 AMPLITUDE PLOT IS PRODUCED. 

IF NOMIT « 1 AMPLITUDE PLOT IS OMITTED. 

INITIAL AMPLITUDES OF F-FUNCTIONS (REMAINING CARDS) I 

AS(«J) IS THE AMPLITUDE OF THE SINE TERM. 

AC(J> IS THE AMPLITUDE OF THE COSINE TERM. 


C 


C 


COMPLEX 

DIMENSION 

1 

2 

3 

4 

5 

6 

7 

8 
9 


YNOZ ( 1 0 > * B(IO)* Cl* C2* C3* CPHITOO)* CSUM* A 
L(10)* N( 10)* S(10)* NAME(IO)* AS(20)* AC(20)* 

U( 250* 40) * AA( 4)* Y(40)* FZ(4*40)* YP(40)* UZ(40)* 

CP( 3* 20* 20)* FRQK20)* DMP1 ( 20) * UMAXf 500)* UAVG( 100)* 
Z(6)* ANGLE(6)* THETA(6)* CFT(6*20)* YI(20)* 

CFTH(6* 20)* CFZ(6* 20)* PRESS(6)* AXVEL(3)* YR(20)* 
TPLOT(SOO)* YPL0T(6* 500)* DUMMYT(500)* DUMMYY(500)* 
IBUF( 3000)* I TT(4)* I TY 1 ( 7 ) * ITY2(7)* ITY3(7)* 

I TY 4 ( 7 ) * ITY5C6)* TAUCUT(20)* ITY6(8)* 

I TP( 3) * TI TLE( 12)* PRS(SOO)* TI(500)* PMAX(SOO)* 
TIMAX(SOO)* YL0(6)* YHI (6)* YLAB( 6) * ITICY(6) 


COMMON 

1 

2 

COMMON 

COMMON 


HV( 20* 4) * C( 3* 20* 20)* D(20*400)* 

KPMAX( 3* 20 )* I C( 3* 20* 20) * KPQMAX(20)* 
IDP(20*400)* I DQ( 20* 400) 

/BLK2/ M( 10)* NS(10)* SJ( 10)* B 

/BLK3/ NJMAX* NLMAX* GAMMA* C0EF(3*20) 


DATA ITT/ ’DIMENSIONLESS TIME* T’/* 

1 I TY1/’ INJECTOR PRESSURE PERTURBATION* THETA =0’/* 

2 ITY2/’ INJECTOR PRESSURE PERTURBATION* THETA = 45V* 

3 I TY3/’ INJECTOR PRESSURE FERTURBATION* THETA = 90V* 

4 I TY4/ ’NOZZLE PRESSURE PERTURBATION* THETA = OV* 

5 I TY 5/ 'NOZZLE AXIAL VELOCITY* THETA = O'/* 

6 I TY6/ ’NOZZLE B.C. ( RE( -GAMMA*Y*PHI T) ) AT THETA * OV* 
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7 I TP/ ' FRESSURE PEAKS * / 

C 

LAST = 250 
ERR = 0*001 
TDEL = 10*0 
NPT = 0 
AA< 1 ) = 0.0 
AA( 2) = 0.5 
AA( 3) = 0.5 
AAC A) = 1.0 
PI = 3. 14 1 5927 
READ < 5* 5003) NOUTCF 

************* COEFFICIENT INPUT SECTION *******^****************** 

THIS VERSION OF LCYC3D READS THE COEFFICIENT DATA FROM 
A FASTRAND FILE GENERATED BY PROGRAM C0EFFS3D* TO READ 
THIS DATA FROM CARDS* USE READ (5*XXXX) INSTEAD OF 
READ (9* XXXX ) IN THIS SECTION. 

INPUT OF MOTOR PARAMETERS AND NUMBER OF TERMS. 

READ (9* 500 1 ) GAMMA* UE* ZE* ZCOMB* NDROPS* NJMAX 
WRITE (6*6001) GAMMA* UE* ZE* ZCOMB* NJMAX 
IF (NDROPS .EG. 0) WRITE (6*6030) 

IF (NDROPS .EQ. 1) WRITE (6*6031) 

NU * 2 * NJMAX 
JMX = NJMAX/ 2 
RLD = 0.5 * ZE 

WRITE (6*6002) 

INPUT OF DESCRIPTION OF SERIES EXPANSION. 

DO 10 K = 1* JMX 

READ (9*5002) NJ* L(NJ)* M(NJ)* N(NJ)* NS(NJ)* S(NJ)* SJ(NJ)* 

1 NAME(NJ) 

WRITE (6*6003) NAME(NJ)* NJ* L(NJ)* M(NJ)* N(NJ)* NS(NJ)* 

1 S( N J ) * SJ(NJ) 

10 CONTINUE 

WRITE (6*6010) 

DO 15 K = U JMX 

READ (9*5010) J* YNOZ(J)* B(J) i 

WRITE (6*6015) J* YNOZ(J)* B(J) 

NJ = (2 * J) - 1 

YR(NJ) = REAL (YNOZ ( J) ) * 

YI (NJ) = Al MAG (YNOZ ( J ) ) 

YR(NJ+1) = YR(NJ) 

YKNJ+l) = YI(NJ) 

15 CONTINUE 
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C ZERO LINEAR COEFFICIENT ARRAYS. 

DO 20 KC = U 3 
DO 20 NJ = 1* 20 
DO 20 NP = 1* 20 
C(KC*NJ,NP) * 0.0 
CP(KC,NJ,NP) = 0*0 
20 CONTINUE 

ZERO NONLINEAR COEFFICIENT ARRAY. 

DO 30 NJ * 1, 20 
DO 30 NPQ * 1* 400 
D(N JsNPQ) = 0.0 
30 CONTINUE 

INPUT OF LINEAR COEFFICIENTS. 

DO 40 KC « 1> 3 
READ (.9* 5003) KM AX 

IF CNOUTCF .GT. 0) WRITE <6*6004) KC* KM AX 
IF <KMAX .EG. 0) GO TO 40 
DO 45 K = 1* KM AX 

READ (9*5004) NJ* NP* CP(KC*NJ*NP) 

IF CNOUTCF .GT. 0) WRITE (6*6005) KC* NJ* NP* CP(KC*NJ*NP) 
45 CONTINUE 
40 CONTINUE 


INPUT OF NONLINEAR COEFFICIENTS. 

READ (9*5003) NLMAX 

IF (NOUTCF .EQ. 2) WRITE (6*6006) NLMAX 
IF (NLMAX .EQ. 0) GO TO 50 
DO 52 NJ = 1* 20 
KPQMAXCN J) = 0 
52 CONTINUE 

DO 55 K = 1* NLMAX 

READ (9*5005) NJ* NP* NG* DT 

IF (NOUTCF .EG. 2) WRITE (6*6007) NJ* NP* NQ* DT 
KPQMAX(NJ) = KPQMAXCN J) + 1 
KPQ * KPQMAX(NJ) 

IDP(NJ*KPQ) = NP 
I DQ(NJ*KPQ) ■ NQ 
D(NJ*KPQ) * DT 
55 CONTINUE 
50 CONTINUE 

************* PRESSURE COEFFICIENT SECTION *********************** 

CALCULATE SPATIAL COORDINATES FOR PRESSURE COMPUTATION. 

DO 51 NFRES * l, 3 
Z(NPRES) » 0.0 
R THETA * NPRES - l 
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ANGLECNPRES) = R THETA * 45*0 
THETACNPRES) = R THETA * PI/4.0 
ZCNPRES + 3> = ZE 
ANGLECNPRES + 3) = ANGLECNPRES) 

THETACNPRES + 3) = THETACNPRES) 

51 CONTINUE 

CALCULATE COEFFICIENTS FOR PRESSURE TIME HISTORIES. 

DO 53 NPRES = 1* 6 
DO 53 J " 1* JMX 
NP « <2 * J) - 1 
Z1 * ZCNPRES) 

ANG = THETACNPRES) 

CALL PHI CFSC J*Z \» ANG# C1*C2*C3> 

IF CNPRES .EQ. 4) CPHITCJ) = Cl 
CFTCNPRESjNP) = REALCC1) 

CFTCNPRES^NP+l) = -AIMAGCC1) 

CFTHCNPRESjNP) = REAL C C2) 

CFTHCNPRES,NP+1) = -AIMAGCC2) 

CFZ (NPRESjNP) = REAL ( C3) 

CFZ (NPRES..NP+ 1 ) s -AIMAGCC3) 

53 CONTINUE 

OUTPUT OF COEFFICIENTS FOR PRESSURE TIME HISTORIES. 

WRITE C 6* 6020 ) 

DO 56 NPRES = U 6 
WRITE C t* 6014) 

DO 56 J = U NJMAX 

WRITE C 6# 6021 ) ZCNPRES)* ANGLECNPRES) * 

1 CFTCNPRES* J)* CFTHCNPRES# J)* CFZCNPRES*J) 

56 CONTINUE 

************* DATA INPUT SECTION ********************************* 

READ C 5*5000) TITLE 

ZERO INITIAL VALUE AND FREQUENCY ARRAYS. 

5 DO 57 K = 1* NJMAX 
ASCK ) = 0.0 
ACCK) = 0.0 

FRQ1CK) = 0.0 1 

57 CONTINUE 


* 

READ COi'jBUSTI ON AND CONTROL PARAMETERS. 

READ C 5# 5006* END = 300) EN* TAU* H# TSTART* TQUI T 

READ CONTROL NUMBERS. 

READ C 5# 5008) NTEST* JMODE* NLOC* NTERMS* NPZ* NOUT 
JMODE = (2 * JMODE) - 1 
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JPMODE = JMODE ♦ NJMAX 
IF (NOUT .GT. 0) NPT ■ 1 


IF (NOUT .EQ. 0) GO TO 9 
READ DATA FOR SETTING UP PLOTS. 

READ (5*5009) YHI(l), YHI(5)* YLAB(l)* YLAB( 5) 

READ (5*5008) ITICY(l)* ITICY(5)* NFI RST* NOMIT 

************ INITIAL AMPLITUDES SECTION ************************** 

9 DO 58 K = 1* NTERMS 

INPUT INITIAL AMPLITUDES FOR F-FUNCTIONS. 

READ (5*5007) J* AST* ACT 
NJ = (2 * J) - 1 
AS(NJ) = AST 
AC(NJ) * ACT 

CALCULATE FREQUENCY AND DAMPING. 

RL = L( J) 

AX ■ RL * PI/ZE 
AXSQ = AX * AX 
SSQ « S(J) * S(J) 

FRQKNJ) = SQRT( SSQ + AXSQ) 

DMPl(NJ) = 0.0 
FRQKNJ+1) = FRQKNJ) 

DMPKNJ+I) = DMPl(NJ) 

CALCULATE INITIAL AMPLITUDES FOR G-FUNCTIONS. 

IF (FRQl(NJ)) 58* 58* 581 
581 GYRU = GAMMA*YR(NJ)*UE 

GY IF ■ GAMMA*YI (NJ)*FRQ1 (NJ) 

GYRF = GAMMA*YR(NJ)*FRQ1(NJ) 

GYIU * GAMMA*YI(NJ)*UE 
C 

NPRES = 4 

IF (NS( J) .EQ. I) NPRES = 6 
C 

A1 » (1.0+ GYRU)*CFZ(NPRES*NJ+ 1 ) 

1 - GYIF*CFT(NPRES*NJ+ 1) 

A2 = GYRF*CFT( NPRES* NJ+ 1 ) + GYIU*CFZ(NPRES*NJ+ 1) 

A3 ■ -(1.0 + GYRU)*CFZ (NPRES* NJ) ♦ GYIF*CFT(NPRES*NJ) 

A4 « GYRF*CFT(NPRES*NJ) + GYI U*CFZ( NPRES* NJ) 

C 

DET * A1*A1 + A2*A2 

IF (DET .LT. 0.000001) GO TO 583 

R1 - A3*AC(NJ) - A4*AS(NJ) 

R2 ■ >A4*AC(NJ) - A3*ASCNJ) 
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C 

ACCNJ+l) = 
ASCNJ+l) = 
GO TO 58 
583 ACCNJ+l) = 
ASCNJ+l) = 
C 

58 CONTINUE 


(HUAI + R2*A2)/DET 
- ( R2* A 1 - R1*A2)/DET 

-ASCNJ) 

AC (NJ) 


OUTPUT OF INITIAL AMPLITUDES* 

WRITE (6,6016) 

DO 590 J « 1, NJMAX 
IF ( AS( J) 5 59 1, 592, 59 1 

592 IF ( ACC J ) ) 591, 590* 591 

591 WRITE (6,6017) J, DMPKJ), FRQ1 ( J), ACCJ), ASCJ) 

590 CONTINUE 

IF CNTEST • EQ • 0) WRITE (6,6025) 

IF CNTEST .£Q. 1) WRITE (6,6026) 

IF <NPZ • EQ*- 1) WRITE (6,6028) 

IF (NOUT *GE* 1) WRITE (6,6027) 

************* LINEAR COEFFICIENTS SECTION ************************ 

DO 59 KC = 1, 3 
DO 59 NJ « 1, 10 

KPMAX(KC,NJ) = 0 
59 CONTINUE 

IF CNPZ .EQ. 0) GO TO 605 
DO 602 J - 1, JMX 
NJ = (2 * J) - 1 
RL «= L( J) 

AX *= RL * PI/ZE 
AXSQ = AX * AX 
SSQ = S(J) * SCJ) 

OMEGA = SQRTCSSQ + AXSQ) 

TAUCUTCNJ) = 2.0 * PI/OMEGA 
TAUCUTCN J+ 1 ) = TAUCUTCNJ) 

602 CONTINUE 


DO 604 NJ = 1, NJMAX 

DO 604 NP = 1, NJMAX 1 

IF (TAU «GT. TAUCUT(NP) ) CP(3,NJ,NP) * 0-0 

604 CONTINUE 

n 

COMPUTE LINEAR COEFFICIENTS FOR GIVEN VALUES OF EN AND TAU. 

605 DO 60 NJ = 1, NJMAX 
DO 60 NP « 1, NJMAX 
CT « CPC 1,N J,NP) 
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IF CCT) 61, 62, 61 

61 KPMAXC 1,NJ) « KPMAXC i,NJ> ♦ 1 
KP « KPMAXC 1,NJ) 

ICC 1,NJ,KP) = NP 
CC1,NJ,KP> * CT 

62 CT a CPC 2,NU,NP) - EW*CPC 3,NJ,NP> 

IF C CT) 63, 64, 63 

63 KPMAXC 2,NJ) = KFMAXC2,NJ) + 1 
KP « KPMAXC 2,NJ) 

ICC2,NJ,KP) = NP 

CC 2,NJ,KP) = CT 

64 CT a EN * CPC3,NJ,NP) 

IF CCT) 65, 60, 65 

65 KPMAXC 3,NJ) = KPMAXC 3, NJ) + 1 
KP = KPMAXC 3, N«J) 

I CC 3,NJ,KP) = NP 
CC3,NJ,KP) = CT 
60 CONTINUE 

************* STEP- SIZE COMPUTATION ****************************** 

NDIV = 1.0 + TAU/H 
RN a NDIV 
H = TAU/RN 
H6 a H/6.0 

************* initial values section ***************************** 

WRITE C 6, 6008) EN, TAU, GAMMA, UE, RLD 
WRITE C 6, 6009) 

WRITE C6, 6022) CANGLECJ), J a 1,6), CANGLECJ), J = 1,3) 

WRITE C6,6012) 

NPl a NDIV + 1 
DO 70 I * 1, NPl 
NSTEP * I - NPl 
RSTEP = NSTEP 
TIME ■ RSTEP * H 
TICI) a TIME 
DO 75 J * 1, NJMAX 
JP « J .+ NUMAX 
IF CACCJ)) 751, 753, 751 
753 IF CASCJ)) 751, 752, 751 
752 UCI,J) = 0.0 
UCI,JP) = 0.0 
GO TO 75 

751 ARG * FRQICJ) * TIME 
FSIN « SINCARG) 

FCOS * CO SC ARG) 

FEXP a EXPCDMP1CJ)*TIME) 

UCI,J> a C ASC J)*FSIN ♦ ACC J)*FCOS) * FEXP 


I 


U(I.JP) = (IAS(J) * FCOS) - (AC(J) * FSIN)) * FRQKJ) * FEXP 
1 ♦ DMPHJ) * U(I.J) 

75 CONTINUE 

C CALCULATE INITIAL VALUES OF PRESSURE AND VELOCITY. 

DO 704 NPRES = Ii 6 
DO 702 J = 1, NJMAX 
COEF( 1. J) = CFT( NPRES. J) 

C0EF(2. J) = CFTH ( NPRES. J) 

COEF(3.d) = CFZ(NPRES.d) 

702 CONTINUE 

DO 703 J = 1. NU 
Y(J) = U(I.J) 

703 CONTINUE 
UBAR =0.0 

IF < NPRES .GT. 3) UBAR = UE 
UMS = 0.0 

IF ( CNDROPS.EQ. 1 ) .AND. (NPRES»LT«4) ) UMS = UE/(ZE*ZC0MB) 

CALL PRSVELCUBAR. UMS.Y.P. VTH.VZ ) 

PRESS (NPRES)' = P 

IF (NPRES .GT. 3) AXV EL (NPRES - 3) = VZ 

704 CONTINUE 

PRS(I) = PRESS(NLOC) 

C 

C CALCULATE INITIAL VALUES OF NOZZLE B.C. 

CSUM = (0. 0*0.0) 

DO 710 J = 1* JMX 

JP = NJMAX + (2 * J) - I 

FT = Y ( JP) 

GT = Y( JP+ 1 ) 

A = CMPLX(FT.GT) 

CSUM = CSUM + YNOZ(J) * CPHIT(J) * A 
710 CONTINUE 

SUM = REAL (CSUM) 

YPHI = -GAMMA * SUM 

WRITE (6.6011) NSTEP. TIME. (PRESS(J). J = 1.6). 

1 (AXVEL(J), J = 1.3). YPHI 

70 CONTINUE 
C 

WRITE (6.6008) EN. TAU. GAMMA. UE. RLD 

WRITE (6.6022) (ANGLE(J). J = 1.6). (ANGLE(J). J = 1.3) 

C 

c ************* initialize control numbers ************************* 

c 

LINE = 8 
K = 0 
MAXNO = 0 
MAXP = 0 

IF (NOUT »EQ« 0) 60 TO 100 

JPLOT = 0 
TMIN = TSTART 
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TNAX = TSTART + TDEL 
YLO < 1 ) = -YHICl) 

DO 90 J = 2.4 
YHI(J) = YHI(l) 

YL0< J) = YLO ( 1 ) 

YLAB(J) = YLABI 1 ) 

I T I CY < J ) b ITICYC1) 

90 CONTINUE 

YLO< 5) = -YHK5) 

YHI<6) = YHI < 5) 

YLO (6) = YL0<5> 

YLAB(6> » YLABCS) 

ITICY <6> « I T I CY ( 5 > 

************* NUMERICAL CALCULATIONS SECTION ********************* 
100 I = NP1 


C 


RUNGE-KUTTA INTEGRATION SCHEKE. 

105 NSTEP = (I - NP 1 + <LAST - NFl) * K) 
RSTEP * NSTEP 
TINE = RSTEP * H 
Tl(l) = TINE 
BO 110 J = 1, NON AX 
JP = J + NJNAX 
RV ( J# 1 ) = U<I-NDIV*JF) 

RV<J*4) = U< I -NDI V+ 1 j Jp ) 

RV(«J* 2) = 0.375+RVC <J, 1 ) + 0.75*RVCJ,4) 
RV< J# 3) s RV(J^2) 

110 CONTINUE 

DO 120 J = 1* NU 
YCJ> b UCI,J> 

120 CONTINUE 

CALL RHSCNU* UY>YP) 

DO 130 J » 1, NU 
FZ ( 1 > J ) b YP< U) 

130 CONTINUE 

DO 140 II b a»4 
DO 144 J = I,- NU 

UZ(U) = Y<«J) ♦ AA< II) * H * FZCII-1.J) 
144 CONTINUE 


0. 12 5+UCI-NDI V+2>UF) 


CALL RHSCNU* II#UZ>YP) 
DO 148 J b 1 # NU 
FZ(H,J) = YPCJ) 

148 CONTINUE 
140 CONTINUE 

DO 1 SO b | f N U 
U(I+liJ) = Y<J) ♦ CFZC 
150 CONTINUE 


l.»J)+2«0*(FZ< 2# J)+FZ( 3> J) ) 


♦ FZ( 4>«J) ) 


* H6 
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C CALCULATE PRESSURE TIKE HISTORIES* 

DO 154 NPRES = 1# 6 
DO 152 J « 1* NJMAX 
COEF ( 1# J ) = CFT(NPhES*J) 

COEF (2* J ) « CFTH(NPRES* J) 

COEF ( 3* U) = CFZ(NPRES*U) 

152 CONTINUE 
UBAR « 0*0 

IF C NPRES .GT. 3) UBAR = UE 
UMS = 0.0 

IF C ( N DROPS • EQ . 1 ) .AND. (NPRES. LT.4) ) UMS « UE/C Z|.*ZCOMB> 
CALL PRSVEL( UBAR-. UM S* Y* P* VTH* VZ ) 

PRESS(NPRES) = P 

IF (NPRES .GT. 3) AXVEL(NPRES - 3) = VZ 
154 CONTINUE 

PRS(I) « FRESS(NLOC) 

CALCULATE VALUES OF NOZZLE B-C. 

CSUM * (0. 0*0-0) 

DO 650 «J = 1* UMX 
UP « NJMAX + (2 * J) - 1 
FT *= Y ( JP> 

GT « Y(JP+1) 

A = CMFLX(FT*GT) 

CSUM * CSUM -«■ YNOZ(J) * CPHIT(J) * A 
650 CONTINUE 

SUM « REAL (CSUM) 

YFHI « -GAMMA * SUM 


DETERMINE MAXIMA AND MINIMA OF PRINCIPAL MODE- AM* PL I TUDE 
FUNCTION FOR USE IN DETERMINING LIMIT-CYCLE BEHAVIOR. 

IF ( U( I * JPMODE) * U(I+1*JPM0DE)) 170* 170* 160 

170 FDEN « U( I * JPMO DE) - U( 1+ 1* JPMODE) 

IF ( PDEN ) 171* 160* 171 

171 PP = U(I* JPKODE)/PDEN 

PA *= (PP - 1.0) * PP * 0*5 
PB = 1.0 - (PP * PP) 

PC * (PP + 1.0) * PP * 0.5 
KAXNO = MAX NO + 1 

UMAX(MAXNO) « PA*U( I - l* JMODE) ♦ PE*U( I * JMODE) + FC*U( I ♦ 1* JMODE) 
IF (MAXNO .GE. 500) GO TO 250 
160 CONTINUE 

DETERMINE MAXIMUM AND MINIMUM PRESSURE AT LOCATION SPECIFIED 
BY NLOC • 

DPL • PRS( I ) - PRS(I- 1) 

DPS = PRS(I-l) - PRS( I -2) 

IF (DPL*DPS) 173* 173* 175 
173 PNUN = PRSU-2) - PRS(I) 
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2.0*PRS(I-1)) 


PDEN * 2.0 * (PRS(I-2) ♦ PRS< I ) - 
IF (PDEN) 174# 175# 174 

174 PP * PNUK/PDEN 
FA » (PP - 1.0) * PF * 0.5 
PB » 1.0 - (PP * PP) 

PC * (PP ♦ 1.0) * PP * 0.5 
MAXF * MAXP ♦ 1 
PKAX(MAXF) = PA*PRS(I-2) + FB*FRS(I-1) + FC*FRS(I) 

UMAX (KAXF) * TI(I-l) + FP*H 

IF (KAXP *6E. 500) GO TO 250 

175 CONTINUE 

IF (NTEST .EQ. 1) GO TO 155 

IF (TIME .LT. TSTART) GO TO 155 

IF ((NOUT .EO. 0) .OR. (NOUT .GT. 6)) GO TO 156 

************* TIME HISTORY PLOTTING SECTION ********************** 
IF ( TMAX .GT* TOUIT) GO TO 156 

IF ((TIME .GT. TMAX) .OR. (JFLOT «GE. 500)) GO TO 1000 

JPLOT = JFLOT ♦ 1 

FILL TIME ARRAY FOR PLOTTING. 

TPLOT( JFLOT) = TIME 

FILL INJECTOR PRESSURE ARRAYS FOR PLOTTING (THETA = 0# 45# 90) 

DO 1001 J = 1# 3 
YFLOT(J# JFLOT) = PRESS(J) 

1001 CONTINUE 

FILL NOZZLE PRESSURE ARRAY FOR PLOTTING (THETA = 0) 

YPL0T(4# JPLOT) = PRESS(4) 

FILL NOZZLE AXIAL VELOCITY ARRAY FOR PLOTTING (THETA * 0) 

YPL0T(5# JPLOT). = AXVEL(l) 

FILL NOZZLE B.C. ARRAY FOR PLOTTING (THETA ■ 0). 

YFL0T(6# JPLOT) = YPHI 

GO TO 156 

1000 NUM = JPLOT 

PLOT TIME HISTORIES. 

DO 1020 NPLOT = NFIRST# NOUT 
C 

JPLOT * 0 
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C ASSIGN FLOTTING PARAMETERS. 

YMIN » YLOCNFLO T) 

YMAX * YHICNFLOT) 

NTICY * I TI CY(NFLOT) 

DELY » YLAB(NFLOT) 

ELIMINATE POINTS THAT ARE OUT OF THE ORDINATE RANGE. 

DO 1010 J = 1* NUM 

IF CCYFLOTCNPLOT* J) .LT. YMIN) .OR. CYPLGTCNPLQT> J) .GT. YMAX)) 
1 GO TO 1010 
JPLOT = JFLOT + 1 
DUMMYT< JPLOT) = TPLOT(J) 

DUMMYY < JFLO T ) « YPLOTCNFLOT# J) 

1010 CONTINUE 

IF C JPLOT .EQ. 0) GO TO 1020 

GO TO C 101 1* 1012# 1013# 1014* 1015* 1016)# NPLOT 

PLOT INJECTOR PRESSURE AT THETA = 0 DEGREES. 

1011 CALL GRAPHSC 1 BUF* 3000* A# JFLGT* 1 1*NTI CY.TMAX# YMAX* TM IN# YMIN* 

1 ITT* I TY 1*2 1/ 4 1* DUMMYT* DUMMYY* 2.0# DELY* TI TLE) 

GO TO 1020 

PLOT INJECTOR FRESSURE AT THETA = AS DEGREES* 

1012 IF <M< JMODE) .EQ. 0) GO TO 1020 

CALL GRAPHSC I BUF* 3000# A# JFLOT* 1 1*NTI CY.TMAX* YMAX* TMIN* YMIN* 

1 I TT* I TY2# 2 1* A2* DUMMYT* DUMMYY* 2.0# DELY^TI TLE) 

GO TO 1020 

PLOT INJECTOR PRESSURE AT THETA = 90 DEGREES. 

1013 IF (M(JMODE) .EQ. 0) GO TO 1020 

CALL GRAPHSC I BUF* 3000* A* JPLOT* 1 1*NTI CY* TM AX* YMAX* TWIN* YMIN* 

1 I TT* I TY 3* 2 1 * A2* DUMMYT* DUMMYY* 2 . 0* DELY* TI TL E) 

GO TO 1020 

PLOT NOZZLE PRESSURE AT THETA = 0 DEGREES. 

101A CALL GRAPHSC I BUF# 3000* A* JPLOT* 1 1* N TI CY* TM AX* YMAX* TM IN* YMIN* 

1 ITT* I TYA* 21* 39* DUMMYT* DUMMYY* 2.0# DELY* TI TLE) 

GO TO 1020 

PLOT NOZZLE AXIAL VELOCITY AT THETA = 0 DEGREES. 

1015 CALL GRAPHSC IBUF* 3000# A* JFLOT* 1 l.NTICY# TKAX* YMAX* TMIN* YMIN* 

1 ITT* I TY5* 21* 32* DUMMYT* DUMMYY* 2.0# DELY* TI TLE) 

GO TO 1020 

PLOT NOZZLE B.C. AT THETA = 0 DEGREES. 

1016 CALL GRAPHSC I BUF* 3000* A# JFLOT* 1 1*NTI CY* TMAX* YMAX# TMIN* YMIN* 

1 ITT* ITY6* 21* AA# DUMMYT* DUMMYY* 2.0* DELY# TI TLE) 

1020 CONTINUE 
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REASSIGN PLOTTING PARAMETERS FOR NEXT SET OF PLOTS. 

JPLOT = 0 

TMN ■= TMAX 

TMAX = TMAX + TDEL 

************* TIME HISTORY PRINTED OUTPUT SECTION **************** 

156 WRITE <6*601 1) NSTEP* TIME* (FRESSCJ)* J = 1*6)* 

1 < AXV EL ( J ) * «J = 1*3)* YFHI 

LINE = LINE + 1 

157 IF CTIME .GT. TQUIT) GO TO 250 
IF (LINE .LT. 52) GO TO 155 
WRITE (6*6013) 

WRITE (6*6022) (ANGLE(J)* J = 1*6)* (ANGLE(J)* J « 1*3) 

LINE - A 

155 I = I ♦ 1 

IF (I .LT. LAST) GO TO 105 

************* LIMIT-CYCLE SECTION ******************************** 

TEST FOR LIMIT CYCLE. 

K » K ♦ 1 

IF ( (NTEST .EO. 0) .OR. (MAXNO .LT. 80)) GO TO 190 

UTOT « 0.0 

DO 180 J = 0. 3 

UMAX « MAXNO - J 

UTOT = UTOT + ABS( UMAX (UMAX ) ) 

180 CONTINUE 

UAVG(K) * UT0T/4.0 
IF <K .EQ. 1) 60 TO 190 

CHANGE = UAVG(K) - UAVG(K-l) 

ABSCHG = ABS( CHANGE/ UAVG(K)) 

IF (ABSCHG .GT. ERR) GO TO 190 

TV « TIME/2.0 

I TV = TM 

ITM = 2* I TM ♦ 2 

TM » ITM 

TSTART = TM + T START 
TQUIT = TV + TQUIT 
IVIN = TSTART 
TMAX = TSTART + TDEL 
NTEST = 0 

RE-ASSIGN ARRAYS. 

190 DO 200 1=1. NPl 

ILAST = LAST - NPl + I 
PRS(I) = PKS< ILAST) 

TI(I) = TI ( ILAST) 
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oo 200 d « 1# NU 

U(I#J) » U(ILAST#J) 
200 CONTINUE 
60 TO 100 


************* PRESSURE maxima and minima printout **************** 

250 WRITE (6*6023) Z(NLOC)# ANGLE(NLOC)# MAXP 
LINE = 4 

DO 255 JST « 1* MAXP* 8 
JSTART = JST 
JSTOP = JST + 7 

IF (JSTOP .GT. MAXP) JSTOP * MAXP 

WRITE (6*6024) (PMAX(J>* J = JSTART* JSTOP) 

WRITE (6*6024) (TIMAX(d)* J = JSTART* JSTOP) 

WRITE (6*6014) 

LINE = LINE + 3 

IF (LINE .LT. 52) GO TO 255 

LINE = 0 

WRITE (6*6013) 

255 CONTINUE 

IF ((NOUT .EG. 0) .OR. (NOMIT .EG. 1>) GO TO 5 

************* PRESSURE MAXIMA PLOTTING SECTION ******************* 

DETERMINE LARGEST VALUE OF FMAX. 

AMPMAX = 0.0 
DO 260 J = 1# MAXP 

IF (PMAX(J) .LT. AMPMAX) GO TO 260 
AMPMAX ■ FMAX(J) 

260 CONTINUE 

RANGE OF PLOT AND COORDINATE LABELING. 

ITM «= AMPMAX +1.0 
AMPMAX * ITM 

ITM = 1.0 + TIMAX(MAXP)/50.0 
TMAX *= ITM * 50 
DELX ■ TMAX/ 10*0 
DELY = AMPMAX/ 10.0 

ELIMINATE NEGATIVE VALUES. 

JPLOT = 0 

DO 262 d • 1# MAXP 
IF (PMAX(J)) 262# 264# 264 
264 JPLOT » JPLOT ♦ 1 

DUMMY T( JPLOT) « TIMAX(J) 

DUKMYY( JPLOT) « FMAX(J) 

262 CONTINUE 
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C PLOT VALUES. 

CALL GRAPHSdBUF# 3000* 4# JPLOT# 101* 10 1* TMAX# AMFMAX* 0.0* 0. 0* 

1 ITT* I TP* 21* 14* DUKMYT# DUMMYY * DELX* DELY# TITLE) 

GO TO 5 

TURN OFF PLOTTING ROUTINE. 

300 IF (NPT .EO. 1) CALL SHPARG 

************* HEAD FORMAT SPECIFICATIONS ************************* 

5000 FORMAT C12A6) 

5001 FORMAT (4F10. 0*215) 

5002 FORMAT < 51 5* 2F 10. 5* IX* A4) 

5003 FORMAT (15) 

5004 FORMAT C2I5*F15.6) 

5005 FORMAT C3I5*F15.6) 

5006 FORMAT (5F10r0) 

5007 FORMAT (I5#2F10.0) 

5008 FORMAT (715) 

5009 FORMAT (7F10.0) 

5010 FORMAT (I5*4F10»5) 

************* VRI TE FORMAT SPECIFICATIONS ************************ 

6001 FORMAT ( 1HI#9H GAMMA ■ #F5« 3# 5X# 5HUE « #F5.3# 

1 5X# 5HZ E ■= # F8 » 5# 5X# 6HZC0MB » #F5.2# 

2 5X#8HNJMAX « #12//) 

6002 FORMAT (2X#29HNAME J L M N NS# 7X# 3HSMN# 3X# 

1 7HJM( SMN)/) 

6003 FORMAT ( 2X# A4# 51 5# 2F 10. 5) 

6004 FORMAT (1H0#26H NUMBER OF COEFFICIENTS C(# I 1# 10H#NJ#NP) IS# I 5/) 

6005 FORMAT ( 2X# 2HC( # I 1, 1H# , I 2# 1H# # I 2# 4H ) = #F10-5) 

6006 FORMAT (1H0#38H NUMBER OF COEFFICIENTS D(NJ#NP#NQ) IS#I5/) 

6007 FORMAT ( 2X# 2HD( # I 2# 1H# # I 2# 1H# # 1 2# 4H ) ■ #F10*5) 

6008 FORMAT( 1H1#45H COMBUSTION PARAMETERS: INTERACTION INDEX = #F7.5# 

1 12X# 1 IHTIME-LAG = # F7 . 5/2X# 17HM0T0R PARAMETERS: # 19X# 

2 8HGAMMA = #F7.5#23H EXIT MACH NUMBER = #F7.5# 

3 22H ' LENGTH/DIAMETER = #F7.5//) 

6009 FORMAT ( 2X# 18HINI TI AL CONDITIONS//) 

6010 FORMAT ( 1H0# 5X# 1HJ# 8X# 2HYR# 8X# 2HY I # 7X# 3HEPS# 7X# 3HETA//) 

6011 FORMAT ( 2X# I 5# FI 2. 5# 10F10. 5) 

6012 FORMAT ( 1H0) 

6013 FORMAT ( 1H1 ) 

6014 FORMAT ( 1H ) 

6015 FORMAT ( 2X# I 5# 4F10. 5) 

6016 FORMAT (1H1#36H INITIAL CONDITIONS ARE OF THE FORM:// 

1 2X#49HU(I#u) * AC( J)*COS( FREQ* T) + AS( J)* SIN (FREQ* T) ) # 

2 14H * EXP(DAMP*T)///6X# 1HJ#8X# 7HDAMPING# 

3 6X#9HFREQUENCY# 10X# 5HAC(J)# 10X# 5HAS(J)//) 
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6017 FORMAT 

6020 FORMAT 
1 

2 

3 

4 

6021 FORMAT 

6022 FORMAT 
1 

2 

3 

4 

6023 FORMAT 
l 

6024 FORMAT 

6025 FORMAT 

6026 FORMAT 

6027 FORMAT 

6028 FORMAT 

6030 FORMAT 

6031 FORMAT 
END 


<2X#I5#4F15#8/) 

<1H1#46H COEFFICIENTS FOR COMPUTATION OF WALL PRESSURE# 
10H WAVEFORM S///43X# 27HC0EFFI Cl ENTS IN SERIES FORI// 

22X# 5H THETA# 10X#4HTIME# 10X# 5H THETA# 10X# 5HAXIAL/ 

6X# 1HJ# 9X# 1HZ# 3X#9H< DEGREES)# 5X# IOHDERI VATI VE# 
5X# 10HDERIVATIVE# 5X# 10HDERIVATI VE//) 

C2X# I 5# F 10- 3# FI 2. 1# 3F15#7) 

C26X# 17HINJECT0R PRESSURE# 1 4X# 1 5HN0ZZLE PRESSURE# 

12X# 21HN0ZZLE AXIAL VELOCI TY/3X# 4H STEP# 8X# 4H TIME# 

F5#0# 5H DEG##F5#0# 5H DEG##F5*0#5H DEG## 

F5#0#5H DEG##F5#0#5H DEG##F5#0#5H DEG## 

F5#0#5H DEG##F5#0#5H DEG##F5#0#5H DEG## 6X# 4HYPHI //) 

( 1H 1 # 38H PRESSURE MAXIMA AND MINIMA AT: Z = #F5#2# 

11H THETA = # F4# 1/ 19H VALUES COMPUTED: #13//) 

<1H # 7X# 8F1 3#6) 

< 2X//2X# 37HTHE TRANSIENT BEHAVIOR IS CALCULATED#) 

C 2X//2X# 39HTHE LIMIT-CYCLE BEHAVIOR IS CALCULATED#) 

C 2X//2X# 33HTHI S RUN PRODUCES PLOTTED OUTPUT#) 

C 2X//2X# * THE PHANTOM ZONES ARE ELIMINATED# * ) 

( 2X# • DROPLET MOMENTUM SOURCE IS NEGLECTED*/) 

C2X#* DROPLET MOMENTUM SOURCE IS INCLUDED*/) 
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SUBROUTINE PHI CFSCNP* Z* THETA* CT* CTH* CZ ) 

THIS SUBROUTINE COMPUTES THE COEFFICIENTS NEEDED TO 
CALCULATE THE WALL PRESSURE PERTURBATION. 


NP IS THE INDEX OF THE COMPLEX SERIES TERM. 

Z IS THE AXIAL LOCATION. 

THETA IS THE AZIMUTHAL LOCATION. 

CT IS THE COEFFICIENT IN THE SERIES FOR THE TIME 
THE VELOCITY POTENTIAL. 


DERIVATIVE OF 


® I H J;5 THE COEFFICIENT IN THE SERIES FOR THE THETA DERIVATIVE 
OF THE VELOCITY POTENTIAL. 


OF THE™LS2iS I pStI I ^l! HE SERIES F0R ™ E AMAL DERIVATIV * 


COMPLEX Cl* CZ* CAXI * CAXIZ* CRAD* CAZI* CAZITH* 
1 B< 10)* CT* CTH* CZ 

COMMON /BLK8/ MC10)* NSC 10)* SJCIO)* B 

Cl x <0»0* 1.0) 

CZ « CMPLXCZ* 0*0) 

CAXI « C CO SRC Cl * BCNP) * CZ) 

CAXIZ * Cl * BCNP) * CSINHCCI * BCNP) * CZ) 

CRAD * CMPLXC SJCNP)*0«0) 

EM = MCNP) 

AR6 = EM * THETA 
FSIN x SINCARG) 

FCOS * COSCARG) 

AZI = FCOS 

IF CNSCNP ) .EQ. 1) AZI = FSIN 
AZITH = EM * FCOS 

IF CNSCNP) «EQ. 2) AZITH = -EM * FSIN 
CAZI = CMPLXC AZ I* 0. 0) 

CAZITH = CMFLX CAZI TH* 0«0) 


CT = CAZI * CAXI * CRAD 
CTH = CAZITH * CAXI * CRAD 
CZ x CAZI * CAXIZ * CRAD 


RETURN 

END 
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SUBRO UT I N E PRSV EL C UBAR# UM S# Y# P# VTH* VZ ) 

THIS SUBROUTINE COMPUTES THE WALL PRESSURE AND VELOCITY. 

UBAR IS THE LOCAL AXIAL STEADY STATE MACH NUMBER. 

UMS IS THE DERIVATIVE OF THE MACH NUMBER FOR THE CASE 
WHEN DROPLET MOMENTUM SOURCES ARE INCLUDED. 

Y IS THE ARRAY CONTAINING VALUES OF THE MODE- AMFLI TUDE 
FUNCTIONS AND THEIR DERIVATIVES. 

PIS THE VALUE OF THE WALL PRESSURE PERTURBATION. 

VTH IS THE TANGENTIAL COMPONENT OF VELOCITY AT THE WALL. 
VZ IS THE AXIAL COMPONENT OF VELOCITY AT THE WALL. 

DIMENSION YC40)# SUMC4)* SUMSQC3) 

COMMON /BLK3/ NJMAX* NLMAX# GAMMA# C0EF<3#20) 

DO 10 I = 1# 4 
SUMCI) = 0.0 
10 CONTINUE 

DO 20 I = 1# 4 
DO 20 J « 1# NJMAX 
JY «= J 

IF < I .EG. 1) JY » J + NJMAX 
II = I 

IF <1 .EG. 4) II ■ 1 

SUM(I) * SUM< I ) ♦ YCJY) * COEFC 1 1 # J) 

20 CONTINUE 

PLIN = SUMCI ) + UBAR* SUM ( 3) + UMS*SUM<4) 

PNL = 0.0 

IF CNLMAX .EG. 0> GO TO 40 
DO 30 I « 1# 3 
SUMSOCI) « SUM ( I > * SUMCI) 

30 CONTINUE 

PNL * 0.5 * C SUMSQC 2) ♦ SIMSQC3) - SUMSQCD) 

C 

40 P = -GAMMA * CFLIN + PNL) 

VTH = SUMC2) 

VZ « SUM (3) 

C 

RETURN 

END 
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SUBROUTINE RHS(NU, 1 1, U, UP) 


DIMENSION 

COMMON 

1 

2 

COMMON 


U(NU), UF(NU) 

RV(20,4), C(3,20,20), D(20,400), 
KPMAX(3,20), IC(3,20,20), KFGMAX( 20), 

I DP (20, 400)* 1 DQC 20; 400) 

/BLK3/ NJMAX, NLMAX, GAMMA, C0EF(3»20) 


DO 10 NJ > 1, NJMAX 
NJP = NJ ♦ NJMAX 
UP(NJ) = U(NOP) 

SL1 * 0.0 
SL2 = 0.0 
SL3 « 0.0 
SNL ■ 0.0 
MAX « KFMAX( 1,N J) 

IF (MAX .EQ. 0> GO TO 25 
DO 20 KP = 1, MAX 
NP = ICC !,NJ,KF) 

SL1 = SL1 ♦ ( C( 1,NJ,KP) * U(NF) ) 

20 CONTINUE 

25 MAX = KPMAX(2,NJ) 

IF (MAX .EQ. 0) GO TO 35 

DO 30 KP * I, MAX 

NFP * I C( 2,NJ,KP) + NJMAX 

SL2 * SL2 + (C(2,NJ,KF) * U(NPP)) 

30 CONTINUE 

35 MAX = KPMAX ( 3,N J) 

IF (MAX .EQ. 0) GO TO 45 
DO 40 KP = 1, MAX 
NP = I C( 3,NJ,KP) 

SL3 «= SL3 ♦ (C(3,NJ,KP) * RV(NP,II)) 

40 CONTINUE 

45 IF (NLMAX .EQ. 0) GO TO 55 
MAX «= KPQMAX(NJ) 

IF (MAX .EQ. 0) GO TO 55 
DO 50 KPQ = 1, MAX 
NP = I DP (NJ, KPQ) 

NQP = IDQCNJjKPQ) + NJMAX 
SNL » SNL ♦ (D(NJ,KPQ) * U(NP> * U(NQP)) 
50 CONTINUE 

55 UP (NJP) = -( SLI + SL2 + SL3 ♦ SNL) 

10 CONTINUE 
RETURN 
END 
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CGMFILER <FLD=ABS) 

SUBROUTINE GRAPHS! IBUF>NLOC*LDEV»NTOT>NTI CX*NTI CY> 

1 XMAX» YM AX* XMIN»YMIN> ITITLX* I TITLY*LTI TLX*LTITLY*XARRAY* 

2 YARRAY*DELX*DELY* TITLE) 


C- 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


IDENTIFIER 


MEANING 


IBUFJ 

NLOCs 

LDEV: 

NTOT* 

NTICX: 

NTICY: 

xmax: 

YMAX: 

XMIN: 

YKIN: 

ITITLX* 

ITITLY* 


ADDhESS OF BUFFER AREA FOR PLOT OUTPUT 
NUMBER OF LOCATIONS IN BUFFER AREA 0=2000) 
LOGICAL DEVICE NUMBER FOR FLOT 
NUMBER OF POINTS TO BE PLOTTED 
NUMBER OF TIC MARKS ON ABSCISSA 0=2) 

NUMBER OF TIC MARKS ON ORDINATE 0=2) 

UPFER LIMIT OF ABSCISSA DOMAIN 
UFFER LIMIT OF ORDINATE RANGE 
LOVER LIMIT OF ABSCISSA DOMAIN 
LOVER LIMIT OF ORDINATE RANGE 
ABSCISSA LABEL 
ORDINATE LABEL 


TYPE 


INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

INTEGER 

REAL 

REAL 

REAL 

REAL 

FI EL DATA ARRAY 
FI EL DATA ARRAY 


LTITLX* NUMBER OF CHARACTERS IN ITITLX 
LTITLY: NUMBER OF CHARACTERS IN ITITLY 

XARRAYS ABSCISSA POINTS IN TERMS OF XMIN-XMAX COORD’S 
YARRAY: ORDINATE POINTS IN TERMS OF YMIN-YMAX COORD’S 
DELX* INTERVALS OF ABSCISSA TIC MARK LABELING 
IN TERMS OF XMIN-XMAX COORDINATES 
DELY: INTERVALS OF ORDINATE TIC MARK LABELING 

IN TERMS OF YMIN-YMAX COORDINATES 
TITLE* LABEL FOR THE VHOLE RUN 


INTEGER 
INTEGER 
REAL ARRAY 
REAL ARRAY 

REAL 

REAL 

FI EL DATA ARRAY 


DIMENSION IBUF!NLCC)*XARRAY!NTOT)*YARRAY!NTOT)* ITITLX! 1 )* 
1 ITITLY C 1)*YLIT! 100) 

DIMENSION TITLE! I) 


FIXED BASIC PARAMETERS 


LOGICAL ZERO 

DEFINEZER0=NDEC«LT»0« AND.ABS! FPN) «LT • • 5 
1 .OR.NDEC.GT.O.AND.ABS(FPN) .LT.5.* 10«**! “NLEC- 1 ) 
DEFINE DNDEC=NDEC-FLD! 0* 36*ZER0)*NDEC-FLD!0* 36*ZER0) 
DEFINE IFIX(FARG)=INT<FARG+.S) 

DATA J/l/ 

DATA HEIGHT/. 105/ 

DATA INTEC/ l / 

DATA ABSCIS/8./ 

DATA ORDINA/6*/ 

DATA I CODE/- 1/ 
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DATA TOPMAR/1./ 

DATA BOTMAR/ 1.5/ 

REAL LEFMAR 
DATA LEFMAR/ 1.9/ 

DATA RYTMAR/1.1/ 

DATA FACT/1./ 

DATA MAXIS/1/ 

DATA MLINE/1/ 

DATA HTLAB/.105/ 



C 

C 19 INITIAL COMPUTATION OF DERIVED PARAMETERS 

C AND INITIAL PLOTS CALL 

C 20 SKIPS PRELIMINARIES FOR 2ND AND SUBSEQUENT CALLS 
C 



GO TO !19,20),J 
19 YDIT(l) = 3./19. 

TICKLE * HEIGHT/2. 

ROTFAC = - 3*/ 14. * HEIGHT - A./7. * HEIGHT 

STARTL = 6 * HEIGHT ♦ ROTFAC ♦ TICKLE 

SEPLAB * STARTL + 1.5 * HEIGHT 

SYMBLH •= 0.070 

REAL LABSEP 

LABSEP = A. * HEIGHT 

ASTART ■ 2. * HEIGHT 

DO 1 I « 2* 100 

1 YDIT(I) «= YDITCI - 1) + C2 * MOD! 1,2) + 1)/19. 

YDITC100J * YDITOOO) ♦ .5 

CALL PLOTS! IBUF,NLOC,LDEV) 

CALL FACTOR! 1.) 

J *= 2 

CALL SYMBOL !HEIGHT,36 * HEIGHT ♦ 5» 5, HEIGHT, TI TLE, 270*, 72) 
CALL PLOT! 1 ., • . 5# ■ 31 
3 DO 2 I = 1,100 

2 CALL PLOT! 0«.,YDI T! I ) , 3 - M0D!I,2)) 

DO 33 I * 1, 100 

33 YDIT!I) ■ YDI T! I ) - ABSCIS - RYTMAR 

C — — — 

C 

C RESET ORIGIN 

C 



XPAGE = BOTMAR ♦ ORDINA 

GO TO 2019 

20 XPAGE « BOTMAR ♦ ORDINA ♦ TOPMAR 

2019 CALL WHERE! RXPAGE, RYPAGE, FACT) 

YPAGE = RYPAGE - LEFMAR 
CALL PLOT! XPAGE, YPAGE, - 3) 

CALL FACTOR! FACT) 
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c 

c 

C DRAW AXES AND LABELING MAXIS TIMES 

C 

C - 

DO 100 I = 1* MAXIS 
100 CALL MY AXIS 

C - 

C 

C DRAW POINTS* OPTIONAL CENTERLINE* AND PAGE SCI SSORLINE 

C MLINE TIMES 

C 

C 

DO 200 I * 1*MLINE 
200 CALL MYLINE 
RETURN 

C 

C 

C ENTRY POINT SHPARG 

C TERMINATE PLOTTING SEQUENCE 

C 

C 

ENTRY SHPARG 

CALL WHERE< RXFAG E* RYPAGE* I ) 

CALL PLOTC RXPAGE* RYPAGE* 999 ) 

RETURN 

C 

C 

C SUBROUTINE MY AXIS (INTERNAL) 

C 

C — 

SUBROUTINE MY AXIS 

STARTL * 6 * HEIGHT + ROTFAC + TICKLE 
IMAX * IFIXCCYMAX - YMIN ) /DELY) 

TICSEP * 0 RDIN A/ ( ABSCNTI CY ) - 1) 

CALL DEN DEC (YM AX* DELY*NDEC) 

K « 1 

N ■ ( ABSCNTI CY)/IMAX) - 1 ♦ MOD< ABSCNTI CY >* 2) 

DO 9 I * 0* IMAX 
GO TO C 1 1* 12) *K 

11 IFC2 * I *LT • IMAX )G0 TO 12 

CALL AXLABC 0* * I TI TLY*LTI TLY*HTLAB) 

K - 2 

12 FPN « YMAX - I * DELY 

I FCZERO ) FPN * 0. 

TMID * 1. 

XPAGE *= - I * ORDINA/I MAX - .5 * HEIGHT 

IFCFPN) 1 13* 122* 1 18 

113 IFCNDEC - 2)115* 114*112 

114 YPAGE « STARTL 
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GO TO 112 

115 IF<NDEC - 1)117*1 16# 112 

116 YPAGE = STARTL - HEIGHT 
GO TO 112 

117 IF(ABSCFPN) - 100 . > 1 19* 1 16* 1 16 

119 IFCABSCFPN) - 10. ) 120* 121* 121 

120 YPAGE = STARTL - 3 * HEIGHT 
GO TO 112 

121 YPAGE = STARTL - 2 * HEIGHT 
GO TO 112 

122 YPAGE = STARTL - 4 * HEIGHT 
GO TO 112 

118 IFCNDEC - 2)123*116*112 

123 IFCNDEC - 1)125*124*112 

124 IFCFPN - 10. >121* 116# 116 

125 IFCFPN - 10.) 122* 120* 126 

126 IFCFPN - 100. >120* 121*127 

127 IFCFPN- 1000.) 121* 1 16* 128 

128 IFCFPN - 10000. >116* 114* 114 

112 NNDEC » DNDEC 

CALL NUMBER CXPAGE*Y PAGE* HEIGHT* FFN* 270* * NNDEC) 
XPAGE = - I * CORDINA/IMAX) 

DO 10 JJ = 1*N 

YPAGE = TICKLE * TMID 

CALL PLOT CXPAGE*YPAGE* 3) 

YPAGE = YPAGE * C - 1 + I/IMAX * .5) 

CALL PLOTC XPAGE* YPAGE* 2) 

IFC I/IMAX) 110*110*9 
110 YPAGE = 0 

CALL PLOTC XPAGE* YPAGE* 3) 

XPAGE = XPAGE - TICSEP 
CALL PLOTC XPAGE* YPAGE* 2) 

TMID *= .5 
10 CONTINUE 

9 CONTINUE 

K « 1 

IMAX « I FIXC CXMAX - XMIN)/DELX> 

TICSEP * ABSCI S/CNTI CX - 1) 

XPAGE = - ASTART - ORDINA 

CALL DENDECCXMAX* DELX* NDEC) 

DO 28 I = 0* TMAX 

STARTL = - I * ABSCI S/IMAX 

GO TO C 24* 25)*K 

24 IFC 2 * I.LT.IMAX)GO TO 25 

CALL AXLABC 270«* I TI TLX*LT1TLX*HTLAB) 

K = 2 

XPAGE = - ASTART - ORDINA 

25 FPN « XMIN ♦ I * DELX 
IFCZER01FPN = 0. 

IFC FPN >813* 822*818 

813 IFCNDEC - 2)815*81 '1*23 

814 YPAGE = STARTL ♦ 16. /7» * HEIGHT 
GO TO 23 

815 IFCNDEC - 1)817*816*23 
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816 YPAGE = STARTL ♦ 25. /14. * HEIGHT 
GO TO 23 

817 IFtABStFPN) - 100. >8 19* 8 16* 8 16 

819 IFtABStFPN) - 10. >820*821*821 

820 YPAGE = STARTL + 11./ 14. * HEIGHT 
GO TO 23 

821 YPAGE * STARTL ♦ 9./7. * HEIGHT 
GO TO 23 

822 YPAGE = STARTL ♦ 2*/7. * HEIGHT 
GO TO 23 

818 IFtNDEC - 2)823*816*23 

823 IFtNDEC - 1)825*824*23 

824 IFtFPN - 10. >821*816*816 

825 IFCFFN - 10. ) 822*820* 826 

826 IFtFPN - 100. >820*821*827 

827 IFtFPN - 1000- >821*816*828 

828 I Ft FFN - 10000. >8 16*8 14*8 14 

23 NNDEC = DNDEC 

28 CALL NUM BERt XPAGE* YPAGE* HEI GHT* FPN* 270* * NNDEC) 

N = tNTICX/IMAX) - 1 + MODtNTI CX* 2 ) 

DO 26 I = IMAX*0* - 1 
TMID - 1 . 

YPAGE = -If ABSCI S/IMAX 
DO 27 JJ ■ 1*N 

XPAGE = - ORDINA - TICKLE * TMID 

CALL PLOTt XPAGE* YPAGE* 3 ) 

XPAGE = XPAGE + t TICKLE + FLDtO* 36* I *NE*0) * TICKLE) * TMID 
CALL PLOTtXPAGE*YPAGE* 2 ) 

I Ft 1)111* 26* 111 
111 XPAGE = - ORDINA 

CALL PLOTtXPAGE*YPAGE* 3) 

YPAGE « YPAGE + TICSEP 
CALL PLOTtXPAGE*YPAG E* 2) 

TMID * .5 
27 CONTINUE 

26 CONTINUE 
RETURN 

C 

c 

C SUBROUTINE MYLINE t INTERNAL) 

C 



SUBROUTINE MYLINE 

ITOP = IFIXttABSCIS + RYTMAR ♦ .5>/ll. * 99 .) 

IBOT « I FIX t RYTMAR/ 11 . * 99.) 

DO 17 I * 1 * N TO T 

XPAGE « tY ARRAY ( I ) - YMAX)/tYMAX - YMIN) * ORDINA 
YPAGE « tXMIN - XARRAY 1 1 ) ) / 1 XM AX - XMIN) * ABSCI S 
17 CALL SYMBOLtXPAGE* YPAGE* SYMBLH* INTEQ* 270.* I CODE) 
IFtNTICY.GE.O)GO TO 22 
XPAGE « - ORDINA/2 < 

YPAGE ■ - ABSCI S 

CALL PLOTtXPAGE* YPAGE* 3 ) 

DO 18 I « IBOT* I TOP 


18 CALL PLOT(XPAGE,YDIT(I )«3 - M0D(I,2)> 

22 XPAGE » TOPMAB 

YPAGE « - ABSCIS - RYTMAR - .5 

CALL PL0T(XPAGE,YPAGE,3) 

DO 21 I = 1,100 

21 CALL PLOT(XPAGE,YDI T( I ), 3 - M0D(I,2)) 

RETURN 



C 

C SUBROUTINE AXLAB (INTERNAL) 

C 

C 

SUBROUTINE AXLAB( ANGLE, I BCD,NCHARX, HEIGHT) 

DIMENSION IBCD< 7 ) 

LOGICAL S 
INTEGER QSQ/' S’/ 

K « 2 

NCHAR = NCHARX 
S a .FALSE. 

IFCABSC ANGLE) . GT.. 1)G0 TO 30 

XPAGE = - ORDINA/2. - NCHAR * HEIGHT/2 

YPAGE = SEPLAB 

GO TO 31 

30 XPAGE « - ORDINA - LABSEP 

YPAGE » - ABSCIS/2. ♦ NCHAR * HEIGHT/2 

31 L START « 6 * MODCNCHAR, 6 ) - 12 
IFCLSTART.EQ. - 12)LSTART a 24 
LOOK ® NCHAR/ 6 + 1.1 
IF(LSTART.EQ. - 6)G0 TO 13 

I F(FLD( 0, 1 2, ' , S' ) • EQ»FLD(LSTART, 12, IBCD(LOOK) ) )G0 TO 15 
GO TO 14 

13 IF(FLD<0,6, *,* ).NE.FLD(30,6,IBCD(L00K - 1)))G0 TO 14 
I F( FLD< 0,6, * S* ) *NE.FLD(0, 6, IBCD(LOOK) ) )G0 TO 14 

15 NCHAR a NCHAR - 1 
S = .TRUE. 

14 CALL SYM BOL( XPAGE, YPAGE, HEI GHT, I BCD, ANGLE, NCHAR) 

I F ( S ) CALL SYMBOL (999 999 «, 2 * HEI GHT/3, GSQ, ANGLE, 2) 
RETURN 



c 

C SUBROUTINE DENDEC (INTERNAL) 

C 



SUBROUTINE DENDEC ( QMAX, DELO,NDEC) 

IF(INT(ABS(QMAX)).GE. 10)G0 TO 5 
IF(AMOD(ABS(QMAX - DELQ), • 1)*GE..01)G0 TO 7 
NDEC = 1 
RETURN 

5 NDEC = - 1 

RETURN 

7 NDEC = 2 

RETURN 
END 
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APPENDIX E 


USER'S MANUAL FOR THE LINEAR STABILITY 
PROGRAMS: LINSOL AND LSTB3D 


General Description 

Two auxiliary programs, LINSOL and LSTB3D, calculate the linear stability 
characteristics of a cylindrical combustion chamber with distributed combus- 
tion and a conventional nozzle. For given values of the operating parameters 
(i.e., n, xj y> u , and L/d) and a given nozzle admittance (i.e., A and cp) , 
Program LINSOL calculates the growth rate, A, and the frequency, (I), of a 
given acoustic mode. For given values of f Program LSTB3D calculates the 
corresponding values of n and co for neutral stability (a = 0) . These progr am s 
are based on an analytical solution of the linearized version of Eqs. (12). 
After a discussion of the linear analysis, Programs LINSOL and LSTB3D will be 
described. 

Linear Analysis 

For a single acoustic mode, dropping the nonlinear terms in Eqs. (12) 
yields the following linear equation: 



+ C A + 



, d [ A(t ‘ 

nC 3 dt + nC 3 dt J = 


(E-l) 


where A(t) is the unknown complex amplitude function for the mode under con- 
sideration and the coefficients are obtained from Eqs. (C-l) through (C-4) 
by dividing by Cq. Thus the coefficients are complex numbers given by: 



mn 


z 

0 

Z'(z e )Z*(z e ) - J Z"Z*dz 


p z e * 

ZZ dz 

J o 


(E-2) 
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^ 6 
2J u(z)Z'Z*dz + yf ||zz*dz + yYZ(z e )Z*(z g ) 

f» z e * 

I ZZ dz 


'e * 
ZZ dz 


rj sK 3 * 


J e * 

ZZ dz 


where the droplet momentum source has been neglected. When the droplet 
momentum source is included, the y in the second term of Eq. (E-3) is re- 
placed by y + 1 ( see Appendix A) . 

The linear solutions are determined by substituting a solution of the 

form: 


A(t) = 


Ah+±v)t 


(E-5) 


into Eq. (E-l) and separating real and imaginary parts to obtain: 


- C lr + A + ( C 2 r “ n( -'3^A - C p n -0) + C^ne A (AcoS(jjx + ^sin^) (E-6) 


2i w 3 


r CL . + (0 o - nC 0 )o) + nC 0 e“ Ax u)Coso)T 

= _ / l± 3' _3 \ 

^ . r\ ^ -At « 4 ' 


20) + n - nC-,e _AT sinU)T 

2x 3 


(E-7) 


where + iC.^, C 2 = + iC 2i , and is always real. The above 

equations are solved numerically by Program LINSOL to obtain the growth rate, 
A, and the frequency, uj, for given values of n and x. 

The equations describing the neutral stability limits are obtained by 
substituting A = 0 into Eqs. (E-6) and (E-7). Solving the resulting equations 
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n 


(e-8) 


for n and to gives: 


Co + C T / 

_ 2r li/oj 

n “ c^(l - cos out) 

2 - 

w = C + co(nC sin out 
lr 3 


c 2i ) 


which are solved numerically by Program LSTB3D. 


(E-9) 


Program LIMSOL 

Program Structure . A flow chart for Program LINSOL is given in Fig. (E-l) . 
This program consists of the following major sections: (l) input, (2) calcula- 

tion of the coefficients C^, C^, and C^, ( 3 ) iterative solution for A and u>, 
and (4) output. 

Input . The input data required by Program LINSOL includes: (l) a title 

for the run, ( 2 ) the chamber parameters Y> u > L/d, and z /z , ( 3 ) several 
control numbers, ( 4 ) the nozzle admittance, ( 5 ) the mode under consideration, 
and (6) the values of n and T for the cases to be run. This data is described 
in the following table where the location number refers to the columns of the 
card and the following three formats are used: alphanumeric characters (A) , 

integers (i) , and numbers with a decimal point (F) . For the "i" formats the 
values are placed in fields of five locations, while a field of ten locations 
is used with the "f" formats. In either case the numbers must be placed in 
the rightmost locations of the allocated field. 

No. of 


Cards 

Location 


Input Item 

Comments 

1 

1-72 

A 

TITLE 

Title of run. 

1 

1-10 

F 

GAMMA 

Specific heat ratio, y. 


n-20 

F 

UE 

Steady state Mach number 

at nozzle entrance, u . 

e 


21-30 

F 

RLD 

Length-to-diameter ratio 
L/D = z /2. 
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No. of 


Cards 

Location 

Tffie 

Input Item 

Comments 


31-40 

F 

ZC0MB 

Length of combustion zone* 

z /z . 
c / e 


41-45 

I 

NDR0PS 

If 0: droplet momentum 
source neglected. 

If 1: droplet momentum 
source included. 


46-50 

I 

NOZZLE 

If 0: quasi-steady nozzle. 
If 1: conventional nozzle. 

If NOZZLE 

51-55 
= 1 : 

I 

NOPT 

If 1: all coefficients in- 
cluded. 

If 2: imaginary parts 
neglected . 

1 

1-10 

F 

YAMPL 

Amplitude factor of nozzle 
admittance, A. 

11-20 

End of input for NOZZLE 

F 

= 1. 

YPHASE 

Phase of nozzle admittance 

cp. 

1 

1-5 

I 

L 

Axial mode number, l 
(OsLslO). 


6-10 

I 

M 

Tangential mode number, m 
(0 <; M ^ 8). 


11-15 

I 

N 

Radial mode number, n 
(0 <; N ^ 5 ) . 


16-20 

I 

NCASES 

Number of cases to be run 
(NCASES <; 100 ) . 

NCASES 

1-10 

F 

TAU 

Time-lag, t. 


11-20 

F 

EN 

Interaction Index, n. 


The title on the first card should identify the mode under consideration. 
On the second card of input all quantities are the same as those given in the 
input to C0EFFS3D (see Appendix C) except NOPT. NOPT gives the option to 
neglect the imaginary parts of the coefficients and C which are an order 
of magnitude smaller than the corresponding real parts. Neglecting these 
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imaginary parts (NOPT = 2) yields linear solutions consistent with the non- 
linear solutions obtained when the small coefficients are neglected (NEGL = 1 
in input to C0EFFS3D) . The values of n and f for the cases to be run are 
given on a series of NCASES cards. These cards are all read and the values 
of T and n are stored in the arrays TAU(j) and EN(j) before any computations 
are made. 

In addition to the above card input, the acoustic frequencies S are 

mn 

also needed for these calculations. As in Program C0EFFS3D these values are 
given in a DATA statement, which is an integral part of the program. 

Calculation of C 2 , and C^. In this section the coefficients C^, C , 

and appearing in Eqs. (E-6) and (E-7) are calculated using Eqs. (E-2) 
through (E-4). As in Program C0EFFS3D the axial acoustic eigenvalues necess- 
ary for these computations are calculated by Subroutines EIGVAL and FCNS, and 
the integrals of the products of two axial eigenfunctions appearing in Eqs. 
(E-2) through (E-4) are computed by Subroutines AXIAL1 and UBAR. Listings of 
these subroutines are given in Appendix C. 

Iterative Solution for a and U) . Equations (E-6) and (E-7) are of the 

form: 


«) 2 = c lr + f(A><«) 
A = g(A,<o) 


(E-10) 


where the quantity f(AjU)) is small compared to C lr and a is small in most 
cases. Starting with an initial guess of 


(E-ll) 
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Eqs. (E-10) are solved iteratively using the following recursion formulas: 


Vl = C lr + f(A k>V 


A k+1 g ^ A k ,UJ 0 


(E-12) 


At each step of the iteration the quantities AA and p are calculated, where 


AA = | 

A k+i 

- \ 

CP = 1 

\+l 

\ 


(E-13) 


and the computations are terminated when k = 40 or when AA and A^ are lass 

-6 

than e = 10 ~ . The process usually converges in less than 15 iterations. 

Output . The output generated by Program LINSOL consists of a restate- 
ment of the input data followed by the calculated results in tabular form. 

For each case the tabulated results give the values of T and n (TAU and EN), 
the corresponding values of the growth rate A and the frequency tu (LAMBDA 
and OMEGA), and the number of iterations (ITER). When ITER is 40 the last 
values of A and <*> are given followed by the warning message "FAILED TO CONVERGE." 

Sample Input and Output . A sample input for the IT mode is given in 
Table E-l followed by the resulting output in Table E-2. 

Program LSTB3D 

Program Structure . A flow chart for Program LSTB3D is given in Figure 
(E-2). This program consists of the following major sections: (l) input, 

(2) calculation of the coefficients C^, C^, and C , ( 3 ) iterative solution 
for n and «) for neutral stability, and (4) output. 
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Table E-l. Sample Input for LINSOL. 



Table E-2. Sample Output for LINSOL. 

IT MODE. 

DROPLET MOMENTUM SOURCE NEGLECTED 

GAMMA. 1.20 UE . .20 L/D - .50000 ZCOMB « 1.00 

AMPL . .02000 PHASE « 45.0 


TAU 

EN 

LAMBDA 

OMEGA 

ITER 

1.400 

•50000 

-.01789 

1*86593 

7 

1.400 

•58396 

•00000 

1.87005 

7 

1.400 

•60000 

•00339 

1.87078 

7 

1.700 

. 50000 

-.00975 

1.63602 

7 

1.700 

• 54490 

-♦00000 

1*63612 

6 

1.700 

•60000 

•01176 

1.83618 

7 

2.000 

•50000 

-.01537 

1.80691 

6 

2.000 

•57562 

•00000 

1.80410 

8 

2.000 

•60000 

•00487 

1.80322 

8 
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Input . The input data required by Program LSTB3D is basically the same 
as required by Program LINSOL. The first two cards, which give the title of 
the case, the chamber parameters, and the control numbers, are identical in 
content and format to those required by LINSOL. The third card gives the 
mode numbers L, m, and n and is followed by a card giving the nozzle admittance 
if a conventional nozzle is specified. The last card gives the values of t 
for the cases to be run. A detailed description of this input is given below. 

No. of 


Cards 

Location 

Type 

Input Item 

Comments 

1 

1-72 

A 

TITLE 

See input for LINSOL. 

1 

l-4o 

F 

GAMMA, UE, 
RLD, ZCOMB 

See input for LINSOL. 


41-55 

I 

NDROPS, 
NOZZLE, NOPT 

See input for LINSOL. 

1 

1-15 

I 

L, M, N 

See input for LINSOL 

If NOZZLE 

= 1 : 




1 

1-20 

F 

YAMPL, YPHASE 

See input for LINSOL. 

End of input for NOZZLE 

= 1. 



1 

1-10 

F 

TAUMIN 

Smallest value of t. 


11-20 

F 

TAUMAX 

Largest value of t. 


21-30 

F 

DELTAU 

Increment in t. 


The last card gives the values of T which are used in the computation of 
the neutral stability limit. Thus computations are begun for f = TAUMIN, t is 
increased by increments of DELTAU, and computations are terminated when 
T ;> TAUMAX. 

After completion of the computations program control returns to the read 
statement for the nozzle admittance, thus neutral stability curves can be cal- 
culated for several different nozzles for the same set of chamber and mode 
parameters . 

Calculation of C^, C^, and C^. The calculation of the coefficients C^, 
and appearing in Eqs. ( e- 8) and (E~9) is performed in the same manner as 
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given in program LINSOL. 


Iterative Solution for n and 0) . The values of n and cu for neutral stabil- 
ity are calculated for each value of t by solving Eqs. (E-8) and (E-9) using 
the following iteration scheme: 

C 0 + C, ./<*>, 

„ _ 2r lr k 

n k 

C^(l - cosCD^t) 

(E-l4) 

“k + l = C lr + \(\ C 3 si ™k‘ T - C 2i> 

The iteration is started by using = /(^ ’ and is stopped when k = 40 or 
An and fti) are less than e = 10 Convergence is usually obtained in less 
than 20 iterations. 

Output . The output generated by Program LSTB3D consists of a restatement 
of the input data followed by the calculated results in tabular form. For 
each value of T in the range TAUMIN sf s TAUMAX, the tabulated results give 
the value of T (TAU), the corresponding values of n and co for neutral stabil- 
ity (EN and OMEGA), and the number of iterations (ITER). If ITER is 40 the 
last values of n and (c computed are given followed by the warning message 
"FAILED TO CONVERGE." 

Sampl e Input and Output . A sample input for the IT mode is given in 
Table E-3 and is followed by the resulting output in Table E-4. 


Table E-3. Sample Input for LSTB3D. 


JLT M o\s> 

1 J 1 i i « 

7 1 * I* II IJ I) 14 1 

£ . 

t4 i; II 1 

INI 

w >1 » i) » n m v n » m n n u u « w « - 


1* *• 41 

n 

4) 1) 

n 

44 41 44 47 41 

•* M Sl U II li H 14 

1 1 1 < j # 

i-.Z O. 

z 

mi 

» 1) » M » J» » 

±. 

m 

4} 4) 

rr 

44 4* 44 4 

Moll 

*4 4* M II M 1) 14 » * 

rui 1 1 1 iii 

rV,*, « V 

1 J 

' 1 ♦ l« 11 1) 11 14 1) 

LLLL 

U u l| 1 

n " « P »• » 1 * » 1 * *» 11 n u j, „ „ „ m 

!•« II 41 II 44 41 44 1 

1 • 

n 

# « 11 

IT 

M S) 14 

JJT 

» M 

I 

LJ_ _ 

1 1 j * ) i 

|o|. \0 

2 

sj 

— 1 

MI1TTT 

LLLL 

mrrnfr 

* 41 4 

IT 

4) 41 

IT 


' 44 44 10 SI 

Ml 1 

»l II 1 1 

□TET 

S X 

J 


LIT 

m n i r~i r~ rf 1 r i '* i *T* ** * w i — i — rT** , M * w * 41 43 <j « « ** «* 

6 J*. 8 o-ifTT 

M 1 II 

IT 

t) im 1 

Tl 

S S4 

I 
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Table E-4. Sample Output for LSTB3D. 

IT MODE. 

DROPLET MOMENTUM SOURCE NEGLECTED 

GAMMA = 1.20 UE ® .20 RLD « .50000 

AM PL = .02000 PHASE = 45-00 


TAU 

EN 

OMEGA 

ITER 

.60000 

1 • 66353 

2.03102 

6 

•70000 

1.31671 

1.99646 

6 

•60000 

1.08462 

1.96911 

6 

•90000 

.92333 

1*94663 

6 

1.00000 

.80765 

1.92753 

6 

1*10000 

•72330 

1*91089 

6 

1.20000 

•66137 

1.89605 

6 

1*30000 

•61616 

1*88255 

6 

1.40000 

• 58 396 

1.87005 

6 

1*50000 

•56230 

1.85827 

6 

1.60000 

•54961 

1.84702 

5 

1.70000 

•54490 

1.83612 

5 

1*80000 

•54769 

1.82542 

6 

1.90000 

.55785 

1*81479 

7 

2.00000 

•57562 

1*80410 

8 

2. 10000 

•60157 

1.79325 

8 

2.20000 

•63666 

1*78210 

9 

2.30000 

•68221 

1.77055 

10 

2.40000 

•74006 

1.75847 

11 

2.50000 

•81258 

1.74575 

13 

2*60000 

•90278 

1*73224 

14 

2.70000 

1*01446 

1.71783 

17 

2.60000 

1.15226 

1.70240 

21 


ZCOMB 


FORTRAN Listings 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


******************* program linsol ******************************* 

THIS PROGRAM COMPUTES THE DAMPING (LAMBDA) AND FREQUENCY 
(OMEGA) FOR GIVEN VALUES OF THE INTERACTION INDEX <EN) AND 
THE TIME-LAG (TAU). THIS PROGRAM IS BASED ON AN ANALYTICAL 
SOLUTION OF THE COMPLEX DIFFERENTIAL EQUATION. 

THE FOLLOWING INPUTS ABE REQUIRED* 

FIRST CARD* 

THE TITLE OF THE CASE. 

SECOND CARD* 

GAMMA IS THE SPECIFIC HEAT RATIO. 

UE IS THE STEADY STATE MACH NU4BER AT THE NOZZLE ENTRANCE. 

RLD IS THE LENGTH- TO-DIAMETER RATIO. 

ZCOMB IS THE LENGTH OF THE COMBUSTION ZONE# EXPRESSED 
AS A FRACTION OF THE CHAMBER LENGTH. 

NDROPS DETERMINES THE PRESENCE OF DROPLET MOMENTUM SOURCES* 

NDROPS * 0 DROPLET MOMENTUM SOURCE NEGLECTED. 

NDROPS - I DROPLET MOMENTUM SOURCE INCLUDED. 

NOZZLE SPECIFIES THE TYPE OF NOZZLE USED* 

NOZZLE • 0 QUASI -STEADY 
NOZZLE « 1 CONVENTIONAL NOZZLE 
NOPT SPECIFIES THE SOLUTIONS DESIRED. 

NOPT - l COUPLING COEFFICIENTS INCLUDED. 

NOPT « 2 COUPLING COEFFICIENTS NEGLECTED. 

THIRD CARD (FOR CONVENTIONAL NOZZLE ONLY)* 

YAMPL IS THE AMPLITUDE OF THE NOZZLE ADMITTANCE. 

YPHASE IS THE PHASE OF THE NOZZLE ADMITTANCE. 

FOURTH CARD* 

THE MODE IS SPECIFIED BY THE INDICES L, M, AND N. 

L IS THE AXIAL MODE NUMBER AND MUST NOT EXCEED 10. 

m tf St ^™ UTHAL MODE NUMBER AND MUST NOT EXCEED 8* 

N IS THE RADIAL MODE NUMBER AND MUST NOT EXCEED 5. 

NCASES IS THE NUMBER OF CASES TO BE RUN. 

REMAINING CARDS* 

TAU IS THE TIME LAG. 

EN IS THE INTERACTION INDEX. 

****************,,*, w ,,*, w * w ^ w ^ ww ^ w ^^^ w ^ 


COMPLEX 

I 

DIMENSION 


1 

2 

3 


REAL 

COMMON 


B 


YNOZ# RESULT# B(10)# BC# AX(4)« Cl# CZE# 

TITLE(72)' ZEPI# ZEP2 ' CC ' CD ' CE# CSSQ ' CAX 
RJROOTC 10# 5)# 

D(5)« OMEGAOOO)# 

EN(IOO)# TAU(IOO) 

LAMBDAOOO) 
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c 

c ************* DATA INPUT SECTION ********************************* 
C 

ERR * 0*000001 
PI « 3* 141 5927 
Cl » <0. 0/1.0) 

c 

C INPUT ROOTS AND VALUES OF BESSEL FUNCTIONS. 

, DATA ((RJROOTCI/J)/ J * 1/5)/ I « 1/9)/ 

1 3*83171/ 7.01559/ 10.17347/ 13-32369/ 16.47063/ 

2 1.84118/ 5*33144/ 8*53632/ 11*70600/ 14*66359/ 

3 3.05424/ 6.70613/ 9*96947/ 13*17037# 16.34752/ 

4 4.20119/ 8.01524/ 11.34592/ 14*58585/ 17.78875/ 

5 5*31755/ 9.28240/ 12.68191/ 15*96411/ 19* 19603/ 

6 6*41562/ 10*51986/ 13*98719/ 17*31284/ 20*57551/ 

7 7*50127/ 11*73494/ 15.26818/ 18.63744/ 21.93172/ 

8 8 * 57784/ 12*9 3239/ 16 * 529 37/ 19*94185/ 23*268 05/ 

9 9*64742/ 14*11552/ 17*77401/ 21*22906/ 24*58720/ 

C 

C INPUT PARAMETERS. 

READ <5/5000) (TITLECI)/ I * 1/ 72) 

READ (5/5001) GAMMA/ UE/ RLD/ ZCOMB/ NDROPS/ NOZZLE/ NOPT 
IF (NOZZLE .EQ. 1) GO TO 5 
C COMPUTE ADMITTANCE FOR QUASI -STEADY NOZZLE. 

YAMPL * (GAMMA - 1*0) * UE/(2*0 * GAMMA) 

YPHASE « 0.0 
GO TO 7 

5 READ (5/5002) YAMPL/ YPHASE 
7 READ (5/5003) L/ M/ N/ NCASES 
C 

THETA « YPHASE * PI/ 180*0 
YR - YAMPL * COS( THETA) 

YI « Y/*1PL * SIN( THETA) 

YNOZ « CMPLX(YR/YI) 

C 

ZE * 2*0 * RLD 
CZE « CMPLX(ZEzO.O) 

CGAM « CMPLX< GAMMA/ 0*0) 

CAX ■ CGAM 

IF (NDROPS *EQ* 1) CAX « CGAM + <l*0/0*0> 

C 

DO 10 J e 1/ NCASES 

* READ (5/5002) TAU(J)/ EN(U) 

10 CONTINUE 

C 

* C ************* PRELIMINARY CALCULATIONS *************************** 
C 

C 

C ASSIGN ARRAYS FOR ROOTS OF BESSEL FUNCTIONS. 

IF (<M *EQ* 0) .AND* (N *EQ. 0)) GO TO 15 
MM « M ♦ 1 
NN - N 

SMN - RJROOT(MMzNN) 
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c 

c 


60 TO 20 
15 SMN » 0*0 

20 SSQ e SMN * SKN 

CSSG * CMPLX(SSG^O.O) 


C 

C 

C 


C 

C 

C 


C 

c 

c 


acousti c eigenvalues. 

BU) f Se^t' SKN#G " !MA ' ZE ' YAMPL ' YPHASE ' RESULT> 

BC • CON JGC RESULT) 

............. CALCULATE INTEGRALS ......................... 

BO 100 NT « I, A 

*' UE ' it -«OHB.RESULT> 

100 CONTINUE 

............. CALCULATE VALUES AT NOZZLE ENTRANCE 

ZEJ * CCO SHC Cl *BC*CZ E) 

ZEP1 * CCO SHC Cl *BC 1 )*CZE) 

ZEP2 * Cl * BC1) * CSINHC CI*BC 1 )*CZE) 

************* CALCULATE COEFFICIENTS *******************++++++++++ 

CD l CCAxHi?^ 0 " AX<2> * ZEP2*ZEJ)/AXC 1) 

CD C CAX+AXC 3) + C 2 • 0* 0«0)*AXC A) 

1 ♦ CGAM*YN0Z*ZEP1*ZEJ)/AXC1) 

CE « CGAK*AXC3)/AXC1> 


C 

C 

C 


D(l> « REALCCC) 

DC 3) * REAL CCD) 

DCS) * REALCCE) 

IF CNOPT . EQ. 2) GO TO 50 
DC2) * AIKAGCCC) 

DCA) « AIMAGCCD) 

GO TO 55 
50 DC 2) * 0*0 
DCA) « 0*0 


****** CALCULATION OF DAMPING AND FREQUENCY 


55 WRITE <6>6001) 
IF CNDROPS .EQ 
IF CNDROPS .EC 
IF CNOPT .EQ. 
WRITE C6«6002) 
IF C NOZZLE .EQ 
WRITE C6«6005> 
WRITE C6«601 1) 
LINE * 1 A 


CTITLECI), 1*1, 7g) 

• 0) WRITE C6* 6020) 

• 1) WRITE C 6# 6021) 

2) WRITE C 6*60 15) 

GAMMA, UE, ELD* ZCOMB 

• 0) WRITE C6«6012) 
YAMPL, yphase 


****************mm 


% 


t 


1 


1 6b 
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CALCULATE INITIAL GUESSES FOR FREQUENCY. 
RL « L 

AXI « RL * PI/ZE 

AXSQ « AXI * AXI 

SSQ * SMN * SMN 

FRO ■ SQRTC SSQ ♦ AXSQ) 


Jt 


ir 


t 


DO 200 J * I, NCASES 
C 

C2R * DC3) - ENCJ) * DCS) 

C3 = ENCJ) * DC 5) 

C 

LAMBDACI) * 0-0 
0MEGAC1) « FRQ 
C 

K * I 

810 X = LAMBDACK) 

Y * OMEGACK) 

XT ■ X * TAUCJ) 

YT «= Y * TAUCJ) 

EX « EXFC-XT) 

SN * SINCYT) 

CS « COSCYT) 

XSQ *= X * X 

WSQ * DC I) «• XBQ ♦ 02R*X - DC4)*Y 
I + C3*EX* CX*CS ♦ Y*SN) 

A ■ DC 2) ♦ C2R*Y ♦ C3*EX*Y*CS 
BB * 2«0*Y ♦ DC A) - C3*EX*SN 
C 

0MEGACK+ 1 ) « SQRTC WSQ) 

LAMBDACK+1) = -A/BB 
C 

IF CK .EQ. 40) GO TO 216 

DX ■ ABSCLAMBDACK+ 1 ) - LAMBDACK)) 

DY * ABSC OMEGACK* 1 ) - OMEGACK ) ) 

K « K ♦ 1 

IF CCDX .LT. ERR) .AND. CDY .LT. ERR)) GO TO 217 
GO TO 210 
C 


216 WRITE C 6# 6009) TAIKJ), 
GO TO 220 
C 


ENCJ), LAMBDACK), OMEGACK), 


217 WRITE C6,6008) TAUCJ), ENCJ), LAMBDACK), OMEGACK), 

220 LINE e LINE ♦ 2 

IF CLINE .LT. 54) GO TO 200 
WRITE C 6, 6007) 

WRITE C6,601 1 ) 

LINE * 4 


200 CONTINUE 


K 


K 
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♦♦*♦*♦**♦**** FORMAT SPECIFICATIONS *♦♦*♦♦*********♦**♦♦********** 


READ FORMATS 

5000 FORMAT (72A1) 

5001 FORMAT (4F10*0*3I5) 

5002 FORMAT <2F10*0) 

5003 FORMAT <41 5) 


VRITE 

6001 FORMAT 

6002 FORMAT 

1 

6005 FORMAT 

6007 FORMAT 

6008 FORMAT 

6009 FORMAT 

6011 FORMAT 
1 

6012 FORMAT 
6015 FORMAT 

6020 FORMAT 

6021 FORMAT 
END 


FORMATS 
< 1H1* 1X*72A1/) 

<2 a 8 8^n? n “ 'F5.2,5 x, 5 HUE = , F5. 2, 5X* 6HL/D - *F8.5* 
5a*6HZC0WB * #F5*2/) 

(lH* j HAMFL “ •» F8 • 5# 5X* 8NPRASE * *F6*1/) 

C 2X/FS. 3*F8» 5* 2F 10 • 5* 1 6/) 

<2X*F5«3*F8»5*2F10»5*I6* 5X* 18HFAILED TO CONVERGE/) 

^ 2X^4111 TER^^^ # 2HEN* AX* 6KLAMBDA* 5X* 5H OMEGA* 

<2X* 19H QUASI -STEADY NOZZLE/) 

(2X* 2AH COUPLING TERMS NEGLECTED/) 

(2X* * DROPLET MOMENTUM SOURCE NEGLECTED*/) 

(2X* * DROPLET MOMENTUM SOURCE INCLUDED*/) 
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************** 41 **** PROGRAM LSTB3D ******************************* 

THIS PROGRAM COMPUTES THE LINEAR STABILITY LIMITS CONSISTENT 
WITH THE THREE-DIMENSIONAL SECOND-ORDER THEORY. 

THE FOLLOWING INPUTS ARE REQUIRED* 

FIRST CARD: 

THE TITLE OF THE CASE. 


SECOND CARD* 

GAMMA IS THE SPECIFIC HEAT RATIO. 

UE IS THE STEADY STATE MACH NUMBER AT THE NOZZLE ENTRANCE. 
RLD IS THE LENGTH- TO-DIAMETER RATIO. 

ZCOMB IS THE LENGTH OF THE COMBUSTION ZONE* EXPRESSED 
AS A FRACTION OF THE CHAMBER LENGTH. 

NDROPS DETERMINES THE PRESENCE OF DROFLET MOMENTUM SOURCES* 
NDROPS - O DROFLET MOMENTUM SOURCE NEGLECTED. 

NDROPS * 1 DROPLET MOMENTUM SOURCE INCLUDED. 

NOZZLE SPECIFIES THE TYPE OF NOZZLE USED* 

NOZZLE = 0 QUASI -STEADY 

NOZZLE * 1 CONVENTIONAL NOZZLE 

NOPT SPECIFIES WHICH SOLUTION WILL BE COMPUTED. 

NOPT « I COUPLING COEFFICIENTS INCLUDED. 

NOPT « 2 COUPLING COEFFICIENTS NEGLECTED. 

THIRD CARD* 

THE MODE IS SPECIFIED BY THE INDICES L. M. AND N. 

L IS THE AXIAL MODE NUMBER AND MUST NOT EXCEED 10. 

MIS THE AZIMUTHAL MODE NUMBER AND MUST NOT EXCEED 8. 

N IS THE RADIAL MODE NUMBER AND MUST NOT EXCEED 5. 

FOURTH CARD <IF CONVENTIONAL NOZZLE): 

YAMPL IS THE AMPLITUDE OF THE NOZZLE ADMITTANCE. 

YPHASE IS THE PHASE OF THE NOZZLE ADMITTANCE. 

REMAINING CARDS* 

TAUMIN IS THE MINIMUM VALUE OF THE TIME-LAG. 

TAUMAX IS THE MAXIMUM VALUE OF THE TIME-LAG. 

DELTAU IS THE INCREMENT IN TIME-LAG. 


I 


C 


****************************************************************** 


COMPLEX 

1 

DIMENSION 

1 

2 

COMMON 


YNOZ, RESULT. B(10>. BC. AX<4). Cl. CZE. 
CGAM. ZEJ. ZEP 1. ZEP2. CC. CD. CE. CSSQ. CAX 
TITLE( 72). 

RJROOTC 10.5). 

OMEGA. 100). EN( 100) 

B 
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************* DATA INPUT SECTION ********************************* 

ERR » 0*000001 

PI - 3* 141 5927 * 

Cl a <0*0* 1*0) 

INPUT ROOTS AND VALUES OF BESSEL FUNCTIONS. 

DATA ((RJR00T(I*U)* J « 1*5)* I = 1*91/ ^ 

1 3*83171* 7*01559* 10*17347* 13*32369* 16*47063* 

2 1*84118* 5*33144* 8*53632* 11*70600* 14*86359* 

3 3*05424* 6*70613* 9*96947* 13*17037* 16*34752* 

4 4*20119* 8*01524* 11*34 592* 14 * 58585* 17*788 7 5* 

5 5*31755* 9*28240* 12*66191* 15*96411* 19*19603* 

6 6*41562* 10*51986* 13*98719* 17*31284* 20*57551* 

7 7*50127* 11*73494* 15*26818* 18*63744* 21*93172* 

8 8*57784* 12*93239* 16*52937* 19*94185* 23*26605* 

9 9*64742* 14*11552* 17*77401* 21*22906* 24*58720/ 

INPUT PARAMETERS. 

READ (5*5000) (TITLE(I)* I » 1* 72) 

READ (5*5001) GAMMA* UE* RLD* ZCOMB* NDROPS* NOZZLE* NOPT 
READ (5*5002) L* M* N 
8 IF (NOZZLE *EQ* 1) GO TO 5 

COMPUTE ADMITTANCE FOR QUASI-STEADY NOZZLE* 

YAMPL * (GAMMA - 1*0) * UE/(2*0 * GAMMA) 

YPHASE *= 0*0 
GO TO 7 

5 READ (5* 5003* END « 300) YAMPL* YPHASE 
7 READ (5*5003* END >300) TAUMIN* T A UMAX* DELTAU 

THETA » YPHASE * PI/ 180*0 
YR « YAMPL * COS( THETA) 

YI - YAMPL * SIN (THETA) 

YNOZ ■ CMPLX(YR*YI ) 

ZE ■ 2*0 * RLD 
CZE « CMPLX(ZE*0*0) 

CGAM b CMPLX( GAMMA* 0*0) 

CAX b CGAM 

IF (NDROPS *EQ* 1) CAX « CGAM * (1. 0*0*0) 

************* PRELIMINARY CALCULATIONS *************************** 

ASSIGN ARRAYS FOR ROOTS OF BESSEL FUNCTIONS. 1 

IF ((M *EQ* 0) .AND* (N *EQ* 0)) GO TO 15 
MM ■ M + 1 
NN - N 

SMN b RJROOT(MM*NN) 

GO TO 20 

15 SMN b 0*0 

20 SSQ « SMN * SMN 

CSSQ - CMPLX( SSQ*0*0) 
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C CALCULATE AXIAL ACOUSTIC EIGENVALUES* 

CALL EIGVAL(L* SMN* GAMMA* ZE*YAMPL* YPHASE* RESULT) 

B< I) = RESULT 
BC « CONJG( RESULT) 

************* CALCULATE AXIAL INTEGRALS ************************** 
DO 100 NT - 1* 4 

CALL AXIAL lCNT* I* 1# UE*ZE*ZCOMB* RESULT) 

AX(NT) « RESULT 
100 CONTINUE 

************* CALCULATE VALUES AT NOZZLE ENTRANCE **************** 

ZEJ - CCOSH(CI*BC*CZE) 

ZEPl ■ CC0SH(CI*B( I)*CZE) 

ZEP2 ■ Cl * B(l) * CSINH(CI*B( 1 )*CZE) 

************* CALCULATE COEFFICIENTS ***************************** 

CC * (CSSQ*AX( I) - AX( 2) + ZEP2*ZEJ>/AX( 1 ) 

CD ■ (CAX*AX( 3) ♦ (2.0»0«0)*AXC4) 

1 ♦ CGAM*YN0Z*ZEP1*ZEU)/AX(1) 

CE - CGAM*AX(3)/AXC 1) 

Cl = REAL ( CC) 

D1 « REAL! CD) 

E « REAL(CE) 

IF <NOPT .EQ. 2) GO TO 50 
C2 > AIMAG(CC) 

D2 * AIM AG CCD) 

GO TO 55 
50 C2 = 0*0 
D2 ■ 0*0 

*********** CALCULATION OF LINEAR STABILITY LIMIT ****** 

55 0MEGAC1) * SQRT(Cl) 

C 

WRITE (6*6001 ) (TITLE(J), J « 1*72) 

IF CNDROPS .EQ. 0) WRITE (6*6025) 

IF (NDROPS .EQ. 1) WRITE (6*6026) 

IF (NOPT .EQ. 2) WRITE (6*6022) 

WRITE (6*6002) GAMMA* UE* RLD* ZCOMB 
IF (NOZZLE .EQ. 0) WRITE (6*6012) 

WRITE (6*6005) YAM PL* YPHASE 
WRITE (6*6010) 

LINE - 12 
C 

TAU * TAUMIN 

370 IF (TAU .GT. TAUMAX) GO TO 8 
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c 

K - 1 

310 VT « OMEGA(K) * TAU 

BB ■ (D1 + C2/0MEGA(K) )/E 

EN<K> = BB/( 1*0 - COS(WT) > X 

G - (E*EN(K)*SIN(VT) - D2) * OMEGA(K) 

OMEGACK+1) - SQRTCC1 ♦ G) 

IF (K .EQ. 40) GO TO 316 ^ 

IF <K .EO. 1) GO TO 311 
DN - ABS(ENCK) - EN(K-l)) 

DW » ABS(OMEGA(K+ 1 ) - OMEGA(K) ) 

IF ((DN .LT. ERR) .AND. <DW .LT. ERR)). GO TO 317 

311 K - K ♦ 1 
GO TO 310 

316 WRITE (6/6013) TAU* ENCK)# OMEGA(K) * K 
GO TO 318 

317 WRITE (6*6014) TAU* EN(K>* OMEGA(K)# K 

318 LINE - LINE ♦ 2 
TAU « TAU ♦ DELTAU 

IF ((LINE .LT. 60) *0R. (TAU *GT. TAUMAX) ) GO TO 370 
WRITE (6*6015) 

WRITE (6*6010) 

LINE « 6 
GO TO 370 

300 CONTINUE 

************* FORMAT SPECIFICATIONS ****************************** 

READ FORMATS 

5000 FORMAT (72A1) 

5001 FORMAT (4F10. 0/315) 

5002 FORMAT (315) 

5003 FORMAT (3F10*0) 

WRITE FORMATS 

6001 FORMAT ( 1H1# IX/ 72A1/) 

6002 FORMAT (2X#8HGAMMA « #F5«2# SX# 5HUE - #F5«2# 5X# 6HRLD » #F8.5# 

1 5X# 8HZC0MB ■ #F5.2/> 

6003 FORMAT ( 2X# A4# 51 5# 4F10* 5/> ) 

6005 FORMAT ( 2X# 7HAMPL « # F8 * 5# 5X# 8H PHASE « ,F7.2/) 

6007 FORMAT ( 1H ) 

6008 FORMAT ( 1H0) X 

6010 FORMAT ( 2X//8X# 3HTAU# 8X# 2HEN# 5X# 5H0MEGA# 6X# 4HI TER/) 

6012 FORMAT (2X# 19HQUASI- STEADY NOZZLE/) 

6013 FORMAT ( 2X# 3F10. 5# 1 10# 5X# 19H FAILED TO CONVERGE/) 

6014 FORMAT (2X#3F10.5#1 10/> 

6015 FORMAT (1H1) 

6022 FORMAT ( 2X# 24HC0UPLING TERMS NEGLECTED/) 

6025 FORMAT (2X# 'DROPLET MOMENTUM SOURCE NEGLECTED'/) 

6026 FORMAT (2X# 'DROPLET MOMENTUM SOURCE INCLUDED'/) 

END 
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